Link-Level Simulator for 5G Localization

We aims to configure all initial parameters for different function modules before executing the simulation. Specifically, it is composed of six types of initial parameters: system parameters, carrier parameters, positioning signal parameters, wireless channel parameters, HI parameters, and ABF parameters. In particular, the simulator considers much more specific parameters to support various 5G deployments. Therefore, each type of parameter is integrated into a class form for ease of configuration. This ensures that the parameters are not only assigned with default values but also easily customized according to different simulation requirements. In addition, the value ranges of these parameters are restricted in the class files to circumvent abnormal configurations.

Specifically, system parameters are composed of the system layout and the basic simulation settings. These parameters should be configured firstly and can thus interact with the rest simulation setting. Carrier parameters refer to the time-frequency resource configurations meeting the NR transmission numerologies. Positioning signal parameters involve the selection of the signal type and the configuration of the particular signal resource set. Wireless channel parameters initialize fine-grained channel modeling parameters. Particularly, the complicated array setting is specified to support the mmWave transmission and the cross-polarized cases. It also contains indicators to activate or deactivate additional modeling components. HI parameters are composed of the indicators to activate or deactivate each of the HI functions as well as the corresponding specific parameters of the HI functions. ABF parameters are composed of the beam sweeping settings and the indicators to activate or deactivate different types of beam sweeping.

We aims to simulate the physical-layer baseband transmission of the positioning signals from signal generation to parameter estimation. The physical layer simulation is divided into five function modules, including positioning signal generation, baseband modulation, receive signal generation, baseband demodulation, and channel parameters estimation. In particular, the ABF and HI functions are appropriately merged into the baseband modulation and receive signal generation modules.

Specifically, according to the selected signal type and the determined signal resource set, positioning signal generation first generates the symbols and the indices of the selected resource set and then maps them to the time-frequency resource grid. In the baseband modulation, the generated resource grid is first converted into the baseband waveform by the orthogonal-frequency-division-multiplexing (OFDM) modulation. Then the ABF, HIs, and transmit power functions at the transmitter are combined into the baseband waveform. In the receive signal generation, the received waveform is obtained by filtering the baseband waveform of each time slot with the CIR as the filter coefficient. In addition, the ABF, HIs, and additive Gaussian noise functions at the receiver are reformulated into the received waveform. Baseband demodulation performs the standard OFDM demodulation and obtains the received resource grid. Channel parameters estimation acquires the channel parameters from the received resource grid according to the known symbols and indices of the reference signal. Herein, the channel frequency response (CFR) is obtained based on the least-square method, and the reference signal received power (RSRP) is obtained according to the signal correlation.

Note that the downsampling and the carrier conversions cannot be validly simulated in the mathematical software and are thus ignored. Fortunately, some HI effects at the carrier frequency band, e.g. the PN effect, can be reformulated at the baseband.

We aims to simulate the wireless channels in various 5G scenarios according to the NR channel model stated in 3GPP TR 38.901. The simulation is composed of four modules: system initialization, large-scale parameters generation, small-scale parameters generation, and multipath coefficients generation. In particular, the time-variant spatial consistency is considered in the simulation to approach more realistic channel characteristics. In addition, to meet the different simulation requirements, the channel model is realized in three cases: light-of-sight (LOS) only, drop-based, and segment-based.

In order to facilitate research on the CSI-based estimation method, several specific display functions are developed in the simulator. The display functions consist of system layout, power-delay-profile (PDP) distribution, CIR distribution, CFR pattern, time-frequency resource pattern, cumulative distribution function (CDF) distribution, and so on. In particular, the orientation of the transceiver, the array structure, and the beamforming pattern can be elaborately configured in the system layout function.

MultiBS and multiuser transmission can be executed simultaneously while meeting the NR transmission protocol. For example, considering SRS transmission, multiple users can transmit their own SRS resource sets to the BSs. The BSs distinguish and demodulate different SRSs according to the dedicated coding information such as the SRS scrambling identity and the resource pattern. The number of users supported in this case depends on the available time frequency resource and the configuration of the SRS resource sets.

Both uplink and downlink transmissions are supported in the simulator. When a positioning signal is selected, the transmission direction is then determined by a parameter transformation module. Specifically, SRS determines the uplink transmission, while other signals such as PRS, CSIRS, and SSB determine the downlink transmission. In the wireless channel model, the arrival and departure parameters for uplink and downlink should be swapped when the transmit direction changes.

Since the physical transmission and wireless channel model are in line with NR standards, the carrier frequency of the simulator ranges from 0.5 GHz to 100 GHz. In addition, the large array can be configured for the mmWave transmission. Moreover, the HI effects and the ABF functions have also been considered in the mmWave transmission.

The presented simulator supports standard 5G scenarios such as indoor, indoor factory, urban macro, urban micro, and rural macro. Meanwhile, the scenario-related parameters can refer to 3GPP TR 38.901. In addition, extended scenarios can also be included in the simulator if the necessary channel parameters of a certain scenario are obtained. For example, one can perform channel measurement in a specific area, such as the underground parking lot, and then estimate the channel parameters required for channel modeling. Moreover, spatial consistency is supported for the static and dynamic simulations. In other words, there exists a cross-correlation between the channel coefficients generated from any two adjacent points in a consistent moving track of a user or from two adjacent users.

Overall, with the supported physical-layer transmission models and extended capabilities, the simulator can not only be applied to the angle estimation but also evaluate other CSI-based localization algorithms. These algorithms include but are not limited to fingerprinting, non-LOS (NLOS) identification, resource scheduling, and machine-learning-based methods.

The CSIRS is one of the downlink reference signals in NR systems, which is used for the downlink channel sounding and the beam management. The CSIRS supports up to 32-antenna-port transmission due to its multiplexing method, and thus both the AOA and angle of departure (AOD) can be estimated with the CSIRS. Moreover, to circumvent the conflict with other downlink physical channels or physical signals, the resource allocation of the CSIRS must conform to some specific principles [24].

The SRS is one of the uplink reference signals in NR systems and its main functions are similar to that of the CSIRS. In particular, the SRS supports up to 4-antenna-port transmission due to its special cyclic shift coding, and thus the AOD estimation with the SRS performs limitedly. In addition, the SRS sequence is generated based on the Zadoff-Chu sequence, which ensures a lower peak-to-average power ratio for the SRS transmission.

For positioning enhancement in NR systems, the PRS and its measuring and reporting mechanism are introduced in release-16 of 3GPP TS 38.211. In order to support measuring the TOAs from multiple BSs, the resource muting function has been included to circumvent the conflict of the PRS resource sets among different BSs. However, the PRS only supports single-antenna port transmission, and thus the AOD information can not be estimated with the PRS.

The SSB is designed for the primary time-frequency domain synchronization during the cell search phase in NR systems. Due to its determined pattern in the resource grid, the SSB cannot be used for channel estimation. Although the SSB can be configured for the initial beam establishment in the mmWave transmission, the beamwidth of the sweeping beams with the SSB is usually too large to estimate the signal orientation directly. Instead, more accurate angles can be calculated by beamforming-based angle estimation techniques.

There are two frequency ranges for the carrier in NR systems, i.e. sub-6 GHz and mmWave. In the sub-6 GHz, the estimated CSI contains the AOA or AOD information, and thus one can configure the CSIRS or SRS resource set for the CSI-based angle estimation. In the mmWave transmission, the ABF-based method is the primary angle estimation method, and can be executed with the CSIRS, SRS, PRS, and SSB resource sets.

Furthermore, all the positioning signals are modeled with the different class files in the simulator. For each signal generation, the symbols and indices of the selected signal in the resource grid are obtained according to the specification in 3GPP TS 38.211.

Two classes of initial channel parameters are considered in the simulation, which are the layout and channel configurations. In particular, we integrate the array configuration into the layout configuration for ease of parameter initialization.

In the array configuration, the initial parameters include array type, numbers of array elements in each row and column, patch antenna type, polarization, array orientation, array element spacing, antenna positions, and so on. When these parameters are determined, the power patterns of all array elements are first determined. Then the power patterns are transformed into the polarized field patterns in the local coordinate system according to the transformation models as stated in clause 7.3.2 of 3GPP TR 38.901 [21]. Finally, the polarized field patterns in the global coordinate system are obtained by the coordinate transformation. Herein, seven classes of antenna power patterns are designed for the simulation, which can efficiently support the mmWave array modeling. Furthermore, we also integrate the three-sector polarized field patterns into the array configuration. When the indicator of the three-sector case is activated, the three-sector antenna positions are first generated by the specific position calculation module. Then the polarized field patterns are obtained according to the above process.

This module first calculates the two-dimensional (2-D) and three-dimensional (3-D) distances from the BSs to the users as well as the corresponding ground reflection distances. In the meanwhile, the outdoor-to-indoor (O2I) and LOS states for each transmission link are assigned by the predefined values or the generations with certain probabilities. If the O2I state for a transmission link is activated, the corresponding indoor distance and indoor penetration loss are then calculated according to the relevant initial parameters. In addition, with the assigned LOS states, all channel model parameters for the determined scenario will be selected. Next, the seven large-scale parameters are generated by the cross-correlation models as stated in the Winner II channel model [27]. These parameters are shadow fading (SF), Ricean K factor (KF), delay spread (DS), azimuth angle spread of departure (ASD), azimuth angle spread of arrival (ASA), elevation angle spread of departure (ESD), and elevation angle spread of arrival (ESA). Moreover, the path loss is then generated.

Herein, to obtain correlated values of these large-scale parameters, random Gaussian-distributed variables for each parameter are first generated on a 2-D position grid, which covers the positions of all the transceivers in the network. Then the variable grid for one large-scale parameter is filtered in two dimensions by an exponential filter with normalized filtering coefficients. The filtering coefficients are generated according to the decorrelation distance of the corresponding large-scale parameter. In the implementation, the 2-D filtering is performed by filtering each of the two dimensions sequentially.

Furthermore, in the cases of spatial consistency, there exists the cross-correlation among the O2I states (and LOS states) for users. It can be resolved by the method of using spatial-consistency random variables. Specifically, four normal-distributed variables are first located in each corner of one correlation grid. Then an interpolation method is used to generate the normal random variables corresponding to the points within the grid. In addition, the uniform distribution can be generated from the phases of the normal random variables. One suitable interpolation function is given by

$\displaystyle\eta_{\alpha,\beta}=$

$\displaystyle\sqrt{1-\frac{\beta}{d_{\rm corr}}}(\sqrt{1-\frac{\alpha}{d_{\rm corr}}}\eta_{0,0}+\sqrt{\frac{\alpha}{d_{\rm corr}}}\eta_{1,0}){+}\sqrt{\frac{\beta}{d_{\rm corr}}}(\sqrt{1{-}\frac{\alpha}{d_{\rm corr}}}\eta_{0,1}{+}\sqrt{\frac{\alpha}{d_{\rm corr}}}\eta_{1,1}),$

(1)

where $\eta_{\alpha,\beta}$ denotes the normal variable at the point $(\alpha,\beta)$, $d_{\rm corr}$ denotes the decorrelation distance, and $\eta_{1,0},\eta_{0,0},\eta_{0,1},\eta_{1,1}$ denote the normal variables located at four corners of the correlation grid.

This module aims to generate the small-scale parameters according to the determined large-scale parameters. Path delays, path powers, arrival $\&$ departure angles, random coupling of rays, and cross-polarization ratios (XPRs) are generated sequentially. In the cases of spatial consistency, the small-scale parameters at $t=0$ are first generated according to the method stated above. Then for $t>0$, new parameters will be generated by the time-variant small-scale parameters generation functions as stated in Procedure A of the spatially-consistent modeling in 3GPP TR 38.901.

Similarly, random initial phases, multipath channel coefficients, as well as LOS-path and ground-reflection channel coefficients are generated sequentially. The pathloss and the shadowing are finally applied to the channel coefficients. Three classes of channel coefficients are developed to meet different requirements, i.e., LOS-only coefficients, static coefficients, and dynamic coefficients. Specifically, LOS-only coefficients only consider the large-scale parameters. This case is mainly used for channel calibration. For the static coefficients, the channel coefficients for different positions of the transceivers are generated independently but the spatial consistency can hold for large-scale parameters. For the dynamic coefficients, the spatial-consistency channel coefficients for the moving track of a user are generated. In addition, the time-varying Doppler shift is implemented. This model can be used to explore the effect of channel coefficients of a moving user on localization algorithms.

To qualify as a 3GPP compatible implementation, the channel model needs to be calibrated. In the channel model calibration specified at clause 7.8 of 3GPP TR 38.901, very detailed simulation assumptions for different calibration cases have been described. With the assumptions, many 3GPP contributors presented their simulation results for each of the metrics. As the simulation results presented by different institutions vary for each of the metrics, the mean value of those presented results is thus used as the benchmark. Therefore, we also verify the channel model by comparing the results from the simulator with the mean of the references.

We present several numerical results of the full calibration and the spatial consistency calibration for the NR channel model as in Figure3 and Figure4, respectively. It can be seen from Figure3 and Figure4 that the wireless channel model presented in this simulator meets the requirements as stated in 3GPP TR 38.901.

The array manifold used for angle estimation is severely distorted due to the array baseline perturbation, the antenna mutual coupling, and the antenna pattern error. Although there have been plenty of works focusing on the compensation methods to the array manifold (see [28] and the references therein), limited results have been obtained in practical applications.

We propose an accurate modeling method by integrating the measured antenna phase offset into the generated CFR. Specifically, the real antenna phase offset is first obtained by measuring the real array in an anechoic chamber and 3-D interpolating. Then letting $\psi_{l,k,\theta,\phi}$ denote the phase offset of the $l$-th receive antenna at the $k$-th subcarrier for the azimuth of angle $\theta$ and the elevation of angle $\phi$, the phase offset sequence for the $p$-th path is given by

$\displaystyle\boldsymbol{\psi}_{p,l}=[{\rm e}^{{\rm j}\psi_{l,\frac{-(K-1)}{2},\theta_{p},\phi_{p}}},...,{\rm e}^{{\rm j}\psi_{l,\frac{(K-1)}{2}-1,\theta_{p},\phi_{p}}}]^{\rm T},$

(2)

where $K$ denotes the length of CFR in frequency domain. Letting $\Delta f,\tau_{p,l}$, and $A_{p,l}$ denote the subcarrier spacing, the delay, and the amplitude of the $p$-th path at the $l$-th antenna, respectively, the $p$-th path CFR is given by

$\displaystyle\boldsymbol{f}_{p,l}=A_{p,l}[{\rm e}^{-{\rm j}2\pi\Delta f\tau_{p,l}\frac{-(K-1)}{2}},...,{\rm e}^{-{\rm j}2\pi\Delta f\tau_{p,l}\frac{(K-1)}{2}}]^{\rm T}.$

(3)

Then the CFR containing the APO at the $l$-th receive antenna is given by

$\displaystyle\boldsymbol{h}_{l}=\sum_{p}\boldsymbol{\psi}_{p,l}\odot\boldsymbol{f}_{p,l}.$

(4)

The proposed modeling method can simulate the actual APO accurately compared to the random generation method. Thus, an effective angle estimation method can be designed while eliminating the APO for a specific antenna array, using the proposed simulator. Some achievements have been published using the proposed modeling method, please see [30, 31] for more information.

Among numerous existing physical-layer simulators for wireless communications, the CIR is usually generated by approximating the real delays with a multiple of the sampling interval. However, this method cannot simulate the energy dispersion effect on the received CIR, which hinders LOS-path detection and TDOA estimation. To resolve this issue, we propose to reformulate the real delays into the CFR of every single path which has been included in equations (3) and (4).

Due to the time synchronization error, the start point of the fast-Fourier-transform window does not correspond to the beginning of each OFDM symbol. This will bring the linear phase offset on the frequency dimension. In addition, there exists a nonlinear offset among receive links due to their different group delays and the clock synchronization error. The adverse phenomenon above is referred to as the TO, which can easily be modeled by multiplying a phase offset to the waveform of each receive link in the frequency domain. Herein, the TO of each link can be randomly generated with the truncated Gaussian distribution.

In the mmWave systems, analog beams are formed by the phased array transceivers [29]. Herein, the phase of each RF link is adjusted by a bit-limited digital phase shifter with a constant phase step. In addition, although some advanced calibration methods have been used to achieve very fine beamsteering resolution and gain-invariant phase tuning, the phase error is still up to several degrees. The amplitude error also exists due to its imperfect compensation method. Thus combining the received signals by the phase and amplitude errors will result in a beamsteering error. To analyze the impact of beamsteering error on angle estimation in mmWave transceivers, we model the effects of phase step, phase error, and amplitude error on beamforming weights. Letting $M$, $\psi$, and $\xi$ denote the bit number of the phase shifter, the phase error, and the amplitude error, respectively, then one beamforming weight is given by

$\displaystyle w=10^{\xi/20}{\rm e}^{{\rm j}(2\pi\frac{i}{2^{M}}+\psi)},\leavevmode\nobreak\ i\in\{0,1,...,2^{M}-1\},$

(5)

where $\xi$ and $\psi$ can be randomly generated with the truncated Gaussian distribution.

The CFO represents the mismatch between the carrier frequency of the received signal and the center frequency of the local oscillator, which incurs a frequency offset of the baseband signals. One can model the CFO by multiplying the received baseband signal with a factor ${\rm e}^{{\rm j}2\pi\epsilon n/K}$, where $\epsilon$ denotes the normalized CFO.

Due to the imperfection of the IQ modulator and/or demodulator at the transceivers of 5G NR systems, the signals in the I and Q branches are not strictly orthogonal. In addition, there exist differences between the IQ signals in transmission lines, analog-to-digital transformation, and low-pass filters, and thus the link gains for the IQ signals are different, which is usually frequency-dependent. The two analog impairments stated above are usually referred to as frequency-independent and frequency-dependent IQ imbalances, respectively.

In the simulator, we use the IQ imbalance models for transceivers as presented in [32]. Specifically, letting $y(n)$ and $x(n)$ denote the complex baseband signals before and after the IQ imbalance distortion, respectively, then IQ imbalance formulation at a transmitter is given by

$\displaystyle y(n)=\mu(n)\otimes x(n)+\nu(n)\otimes x^{*}(n),$

(6)

where $\mu(n)=\left(\frac{\alpha-\beta}{2}\right)g^{\rm I}(n)+\left(\frac{\alpha+\beta}{2}\right)g^{\rm Q}(n),$ and $\nu(n)=\left(\frac{\alpha-\beta}{2}\right)g^{\rm I}(n)-\left(\frac{\alpha+\beta}{2}\right)g^{\rm Q}(n).$ Herein, $\alpha$ and $\beta$ are defined by $\alpha=\cos\psi+{\rm j}\xi\sin\psi$ and $\beta=\xi\cos\psi+{\rm j}\sin\psi$, respectively. In addition, it is assumed that the clocks used by the IQ branches have amplitude mismatch $\xi$ as well as phase mismatch $\psi$, and the discrete-time impulse response of the analog filters in the IQ branches are $g^{\rm I}(n)$ and $g^{\rm Q}(n)$, respectively. If there is no frequency dependency, i.e., assuming that $g^{\rm I}(n)=g^{\rm Q}(n)$, then we have

$\displaystyle y(n)=[\alpha x(n)-\beta x^{*}(n)]\otimes g^{\rm I}(n).$

(7)

Correspondingly, the IQ imbalance formulation at a receiver can be obtained similarly.

The PN is generated from the imperfect of the local oscillator, which is composed of a reference oscillator, a voltage-controlled oscillator, and a phase-locked loop. Usually, the PN incurs the time-variant carrier frequency offset in the radio signal generated from the local oscillator, and the PN effect can be represented by its power spectral density (PSD). In the simulator, we implement a multipole multizero PSD model, which is given by

$\displaystyle S(f)=S(0)\prod_{i}\frac{1+(f/f_{{\rm z},i})^{2}}{1+(f/f_{{\rm p},i})^{2}},$

(8)

where, $f_{{\rm z},i}$ and $f_{{\rm p},i}$ denote the $i$-th zero frequency and the $i$-th pole frequency, respectively. In addition, $S(0)$ denotes the PSD with unit dBc/Hz at $f=0$ Hz.

Generally, the implementation method of the PN model is stated as follows: filtering the Gaussian noise by a low-pass filter, then the PN signal can be obtained by the inverse fast Fourier transform of the filtered signal.

The PAN usually performs severe signal distortion in mmWave transmissions. To explore the effect of PAN on the estimated CSI, we employ a memoryless power amplifier model, in which the signal at the output of the non-linear circuit is given by [33]

$\displaystyle y_{\rm PA}=F_{\rm A}(x){\rm e}^{{\rm j}F_{\rm P}(x)}x/|x|$

(9)

where $x$ denotes the input signal and $F_{\rm A}(x),F_{\rm P}(x)$ are the characteristic functions of PAN. Particularly, the definitions of $F_{\rm A}(x)$ and $F_{\rm P}(x)$ are given in the following.

(i) Rapp amplitude-to-amplitude-conversion model:

$\displaystyle F_{\rm A}(x)=\frac{\eta|x|}{\left(1+\left(\eta|x|/A_{\rm sat}\right)^{2\zeta}\right)^{1/2\zeta}}$

(10)

in root-mean-square (RMS) Volts and

(ii) Modified Rapp amplitude-to-phase conversion model:

$\displaystyle F_{\rm P}(x)=\frac{\alpha|x|^{\gamma_{1}}}{\left(1+\left(|x|/\beta\right)^{\gamma_{2}}\right)}$

(11)

in degrees, where $|x|$ denotes the amplitude of the input signal, $\eta$ denotes the small gain signal, $\zeta$ denotes the smoothness factor, and $A_{\rm sat}$ denotes the high-power amplifier input saturation amplitude. In addition, $\alpha,\beta,\gamma_{1}$, and $\gamma_{2}$ are the fitting coefficients.

This module aims to determine the beam sweeping pattern. The specific beam partition method can be designed by combining the relevant information such as beam number, array size, and sweeping ranges of azimuth and elevation angles. In the simulator, we use the grid-of-beam-based method which partitions the sweeping ranges evenly. Specifically, the 3dB beamwidths for two dimensions are first calculated according to the common approximate formulation for the uniform planar array (UPA) or uniform linear array (ULA), i.e. ${BW}=0.886\lambda/(Nd\cos\theta)$, where $N$ and $d$ denote the antenna number and spacing of an array, respectively. Also, $\lambda,\theta$ denote the wavelength and the angle of arrival/departure of a radio frequency signal, respectively. Then the beam covering interval and all the beam angles can be obtained successively while ensuring that the covering interval is smaller than the 3dB beamwidth.

This module aims to transform the obtained beam angles into the corresponding beamforming weights. Usually, it can be implemented using the steering-vector transformation for a UPA or ULA.

This module aims to determine an optimal beam during the beam sweeping or the beamforming-based angle estimation. For the beam sweeping, one can select the optimal beam based on the maximal RSRP. For the angle estimation, some advanced estimation techniques can be used to obtain more accurate beam angles. As stated in the application cases of the next section, we will present two classes of angle estimation methods, namely differential-beamforming-based and auxiliary beam-pair-based methods. In the physical-layer process of the presented simulator, the transmit and receive beamforming weights are respectively multiplied into the corresponding baseband waveforms.

To verify the simulator performance on the localization, we present an indoor localization simulation²²2One stable version of the presented simulator has been uploaded to the GitHub website; please see https://github.com/Group85GP/Group85GP/tree/5G-positioning-simulator/.. The simulation workflow for the user localization is stated as follows.

1. 1.

   Specify the simulation parameters of the system layout, carrier, SRS, channel, HI, and localization function modules for a user, which are summarized in Table IV.

2. 2.

   Obtain the CFR for each link by the physical-layer transmission process from the SRS generation to the CFR estimation.

3. 3.

   Obtain the user location by the AOA estimation and two-dimensional localization modules.

In Step 1), i.e. parameter initialization, all the BSs are placed with equal spacing and the user is randomly dropped in the system layout, which is in line with the scenario assumption in Table 7.2-2 of 3GPP TR 38.901. In addition, it is assumed that the ULA of each BS is composed of directional antennas.

In Step 2), i.e. physical-layer transmission, when a user is dropped in the layout, the user sends its specific SRS to all the BSs. The SRSs from multiple users can be distinguished by scrambling identity coding or different time-frequency resource patterns. Then each BS estimates CFR from the received signal according to the known SRS configuration information.

In Step 3), i.e. location estimation, $2$ best BSs will be first selected from the $12$ deployed BSs according to a certain selection criterion such as the maximum RSRP. The AOAs will be estimated from the CFRs of the $2$ selected BSs using a spatial process method. Then the user location is obtained by fusing the AOAs from the $2$ selected BSs with a localization algorithm such as the least square method.

Although the CFR contains multipath AOA information, the LOS path is usually much stronger than the other paths. Thus, one can acquire the LOS-path AOA from the CFR using a classical one-dimensional multiple signal classification or digital beamforming method. When the multipath signal strength is high, one can combine spatial processing with delay-domain processing to circumvent the multipath effect.

In the localization, the approximate linear least square is used for the multi-angulation localization. Specifically, It is assumed that $\mathbf{p}=[\alpha,\beta]^{\rm T}$ denotes the 2-D position vector of a user, and ${\boldsymbol{\theta}}$ denotes a column vector, which is composed of the estimated AOAs. Then the observation model for localization is given by ${\boldsymbol{\theta}}=\mathbf{g}(\mathbf{p})+{\mathbf{n}}.$ Herein, ${\mathbf{n}}$ is the noise vector , $\mathbf{g}(\mathbf{p})$ denotes the observation equation in which the $i$-th row is given by ${\rm arctan}\left((\alpha-\alpha_{{\rm BS}_{i}})/(\beta-\beta_{{\rm BS}_{i}})\right)$, and $[\alpha_{{\rm BS}_{i}},\beta_{{\rm BS}_{i}}]^{\rm T}$ denotes the $i$-BS position.

Adopting the first-order Taylor approximation to the observation model, one can obtain the final estimation result as follows:

$\displaystyle{\mathbf{p}}=(\mathbf{B}^{\rm T}\mathbf{B})^{-1}\mathbf{B}^{\rm T}({\boldsymbol{\theta}}-\mathbf{g}(\mathbf{p}_{0}))+\mathbf{p}_{0},$

(12)

where $\mathbf{p}_{0}$ denotes the initial vector of $\mathbf{p}$, the $i$-th row of $\mathbf{B}$ is $[(\beta-\beta_{{\rm BS}_{i}})/\gamma^{2}_{i},(\alpha-\alpha_{{\rm BS}_{i}})/\gamma^{2}_{i}]$, and $\gamma^{2}_{i}=(\beta-\beta_{{\rm BS}_{i}})^{2}+(\alpha-\alpha_{{\rm BS}_{i}})^{2}$.

Furthermore, in order to obtain more accurate results, one can iteratively calculate the equation (12) and substitute $\mathbf{p}_{0}$ with the estimated result in the $(i-1)$ iteration.

First of all, we aim to evaluate the effect of the APO on angle estimation. Figure 5(a) presents the angle estimation error profiles versus different actual angles of incidence for signal-to-noise ratio (SNR) within $\{-5,20\}$ dB. In addition, the actual measured error profile in an anechoic chamber is included for comparison. It can be seen that there exist obvious deviations between the error profile with and without the APO model at the large-angle duration, i.e., near $\pm 60^{\circ}$. In addition, the angle error profiles with the APO model are quite similar for different SNR values.

Moreover, it can be seen that the angle error profiles with the APO model approach the measured angle error profile compared to that without the APO model. It means that the APO model proposed in this simulator can simulate the real APO impact efficiently.

Figure 5(b) performs the estimation CDFs with different HI models and it can be seen that the HI models considered in this simulation significantly deteriorate the estimation precision. Specifically, the simulation considering the TO model achieves the worst results because the inter-channel TO is introduced that can destroy the array manifold destructively. In addition, the estimates considering the APO model first present lower performance than that without HIs and then approaches it. This is because the effect of APO on the angle estimation is not obvious when the angle of incidence is small or at some specific areas of the angle of incidence.

Further, one can also be observed from Figure 5(b) that all the CDFs approach around $80\%$ only. It is because the array of the BS is composed of directional antennas of which the field patterns can only cover limited areas. Therefore, when a user is out of the effective coverage area of the selected BS, the AOA can not be estimated correctly.

In the literature, there are two classes of beamforming-based angle estimations, namely the differential-beamforming-based and auxiliary beam-pair-based methods, which are demonstrated in the following.

Differential-beamforming-based method: In the differential-beamforming-based method, two different beam patterns, namely sum beam pattern and differential beam pattern, are used for mathematical comparison. Herein, the sum beam is the normal beam. The first half of beamforming weights for the differential beam pattern are the same as that for the sum one while the other half of weights are the opposite [35]. The fundamental idea of the differential beamforming is stated as follows.

For ease of expression, it is assumed that a BS equips with $N$-element ULA instead of a UPA while a user equips with a single antenna. In addition, let $\lambda$, $d$, and $\theta$ denote the signal wavelength, the interantenna spacing, and the probing direction, respectively. The beamforming weight vector of the sum beam pattern is given by

$\displaystyle{\bf w}_{\rm sum}=[1,{\rm e}^{{\rm j}\frac{2\pi d}{\lambda}\sin(\theta)},...,{\rm e}^{{\rm j}\frac{2\pi d}{\lambda}(N-1)\sin(\theta)}]=[{\bf v}_{1},{\bf v}_{\rm 2}],$

(13)

Correspondingly, the weight vector for the differential beam pattern is given by

$\displaystyle{\bf w}_{\rm diff}=[{\bf v}_{1},-{\bf v}_{\rm 2}].$

(14)

Next, assume that the actual angle of incidence is denoted by $\theta_{\rm act}$, and then the receive signals by the sum beam and the differential beam are respectively given as

$\displaystyle{\bf y}_{\rm sum}={\bf w}_{\rm sum}{\bf h}x+{\bf n}_{1},$

(15)

and

$\displaystyle{\bf y}_{\rm diff}={\bf w}_{\rm diff}{\bf h}x+{\bf n}_{2},$

(16)

where ${\bf h}=h_{\rm coef}[1,{\rm e}^{-{\rm j}\frac{2\pi d}{\lambda}\sin(\theta_{\rm act})},...,{\rm e}^{-{\rm j}\frac{2\pi d}{\lambda}(N-1)\sin(\theta_{\rm act})}]$, $x$ denotes the receive signal, and $\mathbf{n}$ denotes the noise vector. Herein, $h_{\rm coef}$ denotes the real channel coefficient. Dividing the received signals of the differential beam by the sum beam, one can obtain

$\displaystyle\frac{{\bf y}_{\rm sum}}{{\bf y}_{\rm diff}}\approx\frac{{\bf w}_{\rm sum}{\bf h}}{{\bf w}_{\rm diff}{\bf h}}=\frac{1+{\rm e}^{{\rm j}\alpha}}{1-{\rm e}^{{\rm j}\alpha}},$

(17)

where $\alpha={\rm j}\frac{\pi dN}{\lambda}(\sin(\theta)-\sin(\theta_{\rm act}))$. When the angle deviation $(\theta-\theta_{\rm act})$ is small, this deviation is proportional to $[\sin(\theta)-\sin(\theta_{\rm act})]$. Then from the formulations above, one can obtain that the ratio $\frac{{\bf y}_{\rm sum}}{{\bf y}_{\rm diff}}$ is monotonically increasing (or decreasing) with the angle deviation. When this deviation is covered by both the sum beam and the differential beam, the optimal angle can be achieved.

Auxiliary beam-pair-based method: The auxiliary beam-pair-based method is another way of using the monotonous relation from the ratio of the received signals to the angle deviation. The difference herein is that the spatial frequency deviation is used instead of the angle deviation for ease of derivation. In addition, the ratio between the received signal power by two adjacent beams of the initial beam is used to derive the spatial frequency deviation as presented in [36]. Use the previous assumptions and let $\mu=2\pi d\sin(\theta)/\lambda$ denote the received spatial frequency with the initial beam. Then the spatial frequencies for the two adjacent beams are given by $\mu+\eta$ and $\mu-\eta$, respectively. Herein, we have $\eta=2\pi l/N$ where $l=1,...,N/4$.

Moreover, to circumvent weak received power for one possible auxiliary beam, one can also use a three-beam-based method. Specifically, the received power of two adjacent beams is first compared with each other to reserve the maximal one. Then this reserve beam is combined with the initial beam to derive the spatial frequency deviation by transforming to the two-beam-based problem. The best spatial frequency deviations for two-beam-based and three-beam-based methods are $\pi/N$ and $2\pi/N$, respectively, according to the mathematical derivation.

In this simulation, we aim to evaluate the beamforming-based methods. The simulation workflow for the angle estimation is given in the following.

1. 1.

   Specify the simulation parameters of the system layout, carrier, SRS, channel, and beam sweeping function modules for a transmission link, which are summarized in Table VI.

2. 2.

   Obtain the RSRP for each sweeping beam by the physical-layer transmission process from the SRS generation to the RSRP estimation.

3. 3.

   Obtain the optimal beam direction by an angle estimation module.

In the parameters initialization process, the initial beam orientation at the BS is generated randomly with the given beam orientation ranges. Herein, The initial beam orientation can be acquired from the results of the beam management process as long as this actual orientation is within the 3dB beamwidth of the initial beam. During the physical-layer transmission process, the BS sweeps sum and differential beams or auxiliary beams generated according to the determined initial beam and the angle estimation method. In the angle estimation module, one of the three introduced beamforming-based methods is used to achieve optimal orientation. In addition, the SNR in this simulation is defined by the SNR before the ABF processing.

We first present RMSE profiles for the three-beamforming-based methods in two cases. Specifically, Case (i): the azimuth $\&$ elevation angles of the initial beam are randomly generated from the ranges of $[-51,-39]^{\circ}\leavevmode\nobreak\ \&\leavevmode\nobreak\ [84,96]^{\circ}$ with the truncated normal distribution. Case (ii): the corresponding angles are generated from the ranges of $[-57,-33]^{\circ}\leavevmode\nobreak\ \&\leavevmode\nobreak\ [78,102]^{\circ}$, respectively. Herein, the blue lines with marker ‘$\times$’ represent Case (i), and the others for Case (ii). We observe from Figure 7(a) that as the SNR increases, the RMSEs obtained by the beamforming-based methods decrease. In addition, the estimating results at Case (i) are lower than that at Case (ii). This is because the angle deviation between the initial angles and the actual angles is too large at Case (ii) which deteriorates the estimating performance. Furthermore, for Case (i), the two-beam-based method always performs lower results than the other two methods. Moreover, the sum-and-difference beamforming achieves higher RMSE compared to the other two methods at low SNR but achieves much lower RMSE at high SNR. For Case (ii), the sum-and-difference beamforming can always achieve lower RMSE compared to the other two methods. Therefore, one can conclude that when the initial beam orientation is always within the 3dB beamwidth of the actual beam, the two-beam-based method can be selected to angle estimation. When the range of the initial beam orientation is much large and the SNR is relatively large, the sum-and-difference beamforming can be selected.

Next, we present the CDF results of the three beamforming-based methods for the SNR equaling $-15$ dB and $10$ dB using the first case as stated above. It can be seen from Figure 7(b) that the two-beam-based method can always perform better than other methods for different SNRs. In addition, it can also be observed that the estimating accuracy obtained by the two-beam-based method is within $2$ degree in a percentage of $90\%$ at high SNR setting and with the given assumptions. However, all the estimating methods perform worse accuracy at low SNR settings. Herein, at a low SNR setting, the estimating accuracy is over $6$ degree in a percentage of $90\%$ and over $4$ degree in a percentage of $67\%$, which might be worse than the initial angles obtained by beam sweeping. In the mmWave indoor scenarios, the SNR after the ABF process is usually much higher than $-15$ dB. Therefore, one can conclude that beamforming-based estimation is one of the promising angle estimation methods for angle-based precise localization in mmWave indoor scenarios.

Finally, we analyze the impact of beamsteering error on the ABF-based angle estimation. We use the two-beam-based method only since it has been confirmed to be the most stable method according to the demonstration above. Figure 7(c) presents the CDF results using the two-beam-based method for different beamsteering error cases. Herein, the variables $\Delta_{\Omega}$ and $\Delta_{\psi}$ denote the RMS amplitude and phase errors assigned in the simulation, respectively. In addition, 6-bit digital phase shifters are assumed since it is commonly used at many 5G mmWave chipsets. We observe from Figure 7(c) that the beamsteering error hardly affects the estimating accuracy at different SNR settings. It means that the beamsteering error is not the primary factor that deteriorates the estimating accuracy, although it has been confirmed in [29] that the beamsteering error will incur significant degradation in the sidelobe. Overall, one can say that the high-accuracy estimating results are obtained as long as the actual signal direction of arrival is in the 3dB beamwidth of the determined beam and the SNR is relatively high.