Robust Statistical Beamforming with Multi-Cluster Tracking for Time-Varying Massive MIMO (Extended Version)

TL;DR

This paper introduces a low-complexity, robust joint design for channel estimation, beam tracking, and adaptive beamforming in massive MIMO systems, effectively handling mobility, multi-cluster environments, and frequency selectivity.

Contribution

It proposes a per-cluster channel estimation method combined with a statistical beamforming approach, reducing complexity and improving performance in dynamic massive MIMO scenarios.

Findings

Outperforms conventional methods in dynamic environments

Reduces computational complexity significantly

Effective in multi-cluster and frequency-selective channels

Abstract

In this paper, a joint design of instantaneous channel estimation, beam tracking, and adaptive beamformer construction for a massive multiple-input multiple-output (MIMO) system is proposed. This design focuses on efficiency in terms of performance and computational complexity under the adverse effects of time variation and mobility of sources, the presence of multiuser and multipath components, or simply multi-clusters, and the near-far effect. The design is also suitable for hybrid beamforming and frequency-selective channels. In the proposed system, channel parameters are estimated in time-domain duplex (TDD) uplink mode using a per-cluster approach rather than a joint approach, which significantly reduces the complexity. Per-cluster estimation is possible thanks to the proposed interference-aware statistical beamforming method, namely reduced dimensional Generalized Eigenbeamformer (RD-GEB), which undertakes the computational load of interference mitigation and enables a simpler design for the remaining stages. In addition, the overall design is based on the separation of channel parameters as fast-time and slow-time, leaving only the instantaneous channel estimation and channel matched filtering as fast-time operations, which are handled inside cluster-specific reduced dimensional subspaces. Beam tracking and beamformer construction are held in slow-time rarely, which reduces the time-averaged complexity. Furthermore, beam tracking is performed by leveraging a batch of instantaneous channel estimates, which removes the need for an additional training process. The proposed low-complexity design is shown to outperform the conventional methods.

Tables

Table 1. TABLE I: Nomenclature

Variable	Description
$n$, $p$, $k$, $m$, $u$	Index for symbol, FT-CPI, ST-CPI, cluster, user
$N$, $R$, $D_m$	# Antennas, RFCs (ABF outputs), DBF outputs ($\triangleq$ *)
$P$, $M$, $U$	# FT-CPIs in an ST-CPI, clusters, users
$N_{\text{F}}$, $N_{\text{S}}$	# Symbols in (training, data) mode in an FT-CPI,
$𝑺$, ${𝑾}^{\left(m\right)}$, ${𝑻}^{\left(m\right)}$	ABF, DBF, Total BF matrices
${𝒚}_{n,p}$, ${𝒓}_{n,p}$, ${𝒛}_{n,p}^{\left(m\right)}$	Received Signal (at *)
${𝒉}_p^{\left(m\right)}$, ${\overline{𝒉}}_p^{\left(m\right)}$, ${\widetilde{𝒉}}_p^{(m,m)}$	Channel (at *)
$s_{n,p}^{\left(m\right)}$, $b_{n,p}^{\left(u\right)}$	Symbols from clusters and users
$E^{\left(m\right)}$, $N_0$	Cluster power, noise variance
${\alpha }_{l,p}^{\left(m\right)}$, $𝒂\left(\theta \right)$	Channel complex gain, ULA steering vector
${𝑹}_p^{\left(m\right)}$, ${\overline{𝑹}}_p^{\left(m\right)}$, ${\widetilde{𝑹}}_p^{(m,m)}$	Channel covariance matrix (at *)
${𝚿}_p$, ${\overline{𝚿}}_p$, ${\widetilde{𝚿}}_p^{\left(m\right)}$	Total covariance matrix (at *)
${𝒙}_p^{\left(m\right)}$, ${𝝂}_p^{\left(m\right)}$, $𝑨$	Mobility state vector, innovation, transition matrix
${\theta }_p^{\left(m\right)}$, ${\Delta }_p^{\left(m\right)}$	Mean AoA, AS
${\widehat{\tilde{𝒉}}}_p^{\left(m\right)}$, ${𝝃}_p^{\left(m\right)}$	IEC estimate, estimation error
${𝑹}_k^{(\text{ST},m)}$, ${𝚿}_k^{\left(\text{ST}\right)}$, ${\overline{𝑹}}_k^{(\text{ST},m)}$,	Some ST model variables ( ${(\cdot )}_k^{\left(\text{ST}\right)}$ ) corresponding
${\overline{𝚿}}_k^{\left(\text{ST}\right)}$, ${𝒙}_k^{(\text{ST},m)}$, ${\theta }_k^{(\text{ST},m)}$	to the previously defined variables ( ${(\cdot )}_p$ )

Table 2. TABLE II: FT-CPI and ST-CPI

Number of Symbols in an FT-CPI
	Speed $v$ (m/s)
$f_c/W$	0.1	1	10
30	10e6	1e6	100e3
100	3e6	300e3	30e3
300	1e6	100e3	10e3
1000	300e3	30e3	3e3

Table 3. TABLE III: Average Computational Complexity Per FT-CPI

Task	Technique	# of Multiplications	Dominant Op.
Inst. Eff.	BA-LS	$𝒪\left(M{N}_{\text{F}}{D}_m\right)$	(19)
Channel	Joint LS	$𝒪\left(M^3\right)+𝒪\left(M{N}_{\text{F}}R\right)$	MI* & MM
Estimation	Joint MMSE	$𝒪\left(M^3{R}^3\right)+𝒪\left(M^2{N}_{\text{F}}{R}^3\right)$	MI* & MM
	BA-ML	$𝒪\left(M{N}_{\theta }N{D}_m^2/P\right)$	(30), (35)
Beam		$+𝒪\left(M{D}_m^2+M{N}_{\theta }{D}_m^3/P\right)$
Tracking	SEKF	$𝒪\left(MN{D}_m^2/P\right)$	(30), (36), (41)*
		$+𝒪\left(M{D}_m^2\right)+𝒪\left(M{D}_m^4/P\right)$
	OMP	$𝒪\left(M{N}_{\theta }R{N}_{\text{F}}\right)$	Step 6 in Alg. 1
	FD-GEB	$𝒪\left(M{N}^3/P\right)$	ED
Beamformer	RD-GEB	$𝒪\left(M{R}^3/P\right)$	ED
Construction	FD-MMSE BF	$𝒪\left(M^3+N{M}^2\right)$	MI* & MM
	RD-MMSE BF	$𝒪\left(M^3+R{M}^2\right)$	MI* & MM

Table 4. TABLE IV: Simulation Parameters

Parameter	Description	Value / Details
$N$	Number of antennas	128
$M$	Number of clusters	4
$U$	Number of users	4
$N_{\text{F}}$	FT training sequence length	10
$N_{\text{S}}$	Data mode length in an FT-CPI	990
$T_{\text{F}}$	Duration of an FT-CPI	1e-5 seconds
$P$	Number of FT-CPIs	$\sim$ 1000
$R$	Number of RFCs	16
$D_m$	Number of DBF outputs	3
${\theta }_0^{\left(m\right)}$	Initial AoAs	{${10}^{\circ }$,${20}^{\circ }$,$-{10}^{\circ }$,$-{20}^{\circ }$}
${\Delta }^{\left(m\right)}$	ASs	{$3^{\circ }$,$3^{\circ }$,$3^{\circ }$,$3^{\circ }$}
$E^{\left(m\right)}/N_0$	Cluster SNRs	{10,40,30,30} dB

Equations

$${\mathbf{y}}_{n,p}=\sum _{m=1}^M\sqrt{E^{(m)}}{\mathbf{h}}_p^{(m)}{s}_{n,p}^{(m)}+{\mathbf{\eta }}_{n,p}$$

$$s_{n,p}^{(m)}=b_{n-l_m,p}^{(\mathcal{U}(m))},$$

$${\mathbf{h}}_p^{(m)}=\frac{1}{\sqrt{L}}\sum _{l=1}^L{\alpha }_{l,p}^{(m)}\mathbf{a}({\tilde{\theta }}_{l,p}^{(m)})$$

$${\mathbf{R}}_p^{(m)}=\frac{1}{{\Delta }_p^{(m)}}{\int }_{{\theta }_p^{(m)}-{\Delta }_p^{(m)}/2}^{{\theta }_p^{(m)}+{\Delta }_p^{(m)}/2}\mathbf{a}(\theta ){\mathbf{a}}^H(\theta )d\theta$$

$${\mathbf{\Psi }}_p=\sum _{m=1}^M{E}^{(m)}{\mathbf{R}}_p^{(m)}+N_0{\mathbf{I}}_N.$$

$${\mathbf{x}}_{p+1}^{(m)}=\mathbf{A}{\mathbf{x}}_p^{(m)}+{\mathbf{\nu }}_p^{(m)},\mathbf{A}=\left[\begin{array}{cc}1& T_{\text{F}}\\ 0& 1\end{array}\right],$$

$${\mathbf{r}}_{n,p}\triangleq {\mathbf{S}}^H{\mathbf{y}}_{n,p}=\sum _{m=1}^M\sqrt{E^{(m)}}{\bar{\mathbf{h}}}_p^{(m)}{s}_{n,p}^{(m)}+{\bar{\mathbf{\eta }}}_{n,p}$$

$${\bar{\mathbf{R}}}_p^{(m)}={\mathbf{S}}^H{\mathbf{R}}_p^{(m)}\mathbf{S},$$

$${\bar{\mathbf{\Psi }}}_p={\mathbf{S}}^H{\mathbf{\Psi }}_p\mathbf{S}=\sum _{m=1}^M{E}^{(m)}{\bar{\mathbf{R}}}_p^{(m)}+N_0{\mathbf{I}}_R,$$

$${\mathbf{z}}_{n,p}^{(m)}=({\mathbf{T}}^{(m)}{)}^H{\mathbf{y}}_{n,p}=\sum _{m^{\prime }=1}^M\sqrt{E^{(m^{\prime })}}{\tilde{\mathbf{h}}}_p^{(m,m^{\prime })}{s}_{n,p}^{(m^{\prime })}+{\tilde{\mathbf{\eta }}}_{n,p}^{(m)}$$

$${\tilde{\mathbf{R}}}_p^{(m,m^{\prime })}=({\mathbf{T}}^{(m)}{)}^H{\mathbf{R}}_p^{(m^{\prime })}{\mathbf{T}}^{(m)}$$

$${\tilde{\mathbf{\Psi }}}_p^{(m)}=({\mathbf{T}}^{(m)}{)}^H{\mathbf{\Psi }}_p{\mathbf{T}}^{(m)}=\sum _{m^{\prime }=1}^M{E}^{(m^{\prime })}{\tilde{\mathbf{R}}}_p^{(m,m^{\prime })}+N_0{\mathbf{I}}_{D_m}$$

$${\hat{s}}_{n,p}^{(m)}\triangleq ({\widehat{\tilde{\mathbf{h}}}}_p^{(m,m)}{)}^H{\mathbf{z}}_{n,p}^{(m)}$$

$${\hat{s}}_{n,p}^{(m)}=\sum _{m^{\prime }=1}^M\sqrt{E^{(m^{\prime })}}({\widehat{\tilde{\mathbf{h}}}}_p^{(m,m)}{)}^H{\tilde{\mathbf{h}}}_p^{(m,m^{\prime })}{s}_{n,p}^{(m^{\prime })}+({\widehat{\tilde{\mathbf{h}}}}_p^{(m,m)}{)}^H{\tilde{\mathbf{\eta }}}_{n,p}^{(m)}$$

$${\hat{s}}_{n,p}^{(m)}=$$

$$+({\widehat{\tilde{\mathbf{h}}}}_p^{(m,m)}{)}^H{\tilde{\mathbf{\eta }}}_{n,p}^{(m)}$$

$${\tilde{b}}_{n,p}^{(m)}\triangleq \frac{{\hat{s}}_{n+l_m,p}^{(m)}}{\sqrt{E^{(m)}}({\widehat{\tilde{\mathbf{h}}}}_p^{(m,m)}{)}^H{\widehat{\tilde{\mathbf{h}}}}_p^{(m,m)}}$$

$${\hat{b}}_{n,p}^{(u)}\triangleq \sum _{m=1}^M{p}_m^{(u)}{\tilde{b}}_{n,p}^{(m)}\cong \left(\sum _{m=1}^M{p}_m^{(u)}\right)b_{n,p}^{(u)}+\sum _{m=1}^M{p}_m^{(u)}{e}_{n,p}^{(m)}$$

$${\mathbf{z}}_{:,p}^{(m)}=\sum _{m^{\prime }=1}^M\sqrt{E^{(m^{\prime })}}({\mathbf{s}}^{(m^{\prime })}\otimes {\mathbf{I}}_{D_m}){\tilde{\mathbf{h}}}_p^{(m,m^{\prime })}+{\tilde{\mathbf{\eta }}}_{:,p}^{(m)},$$

$${\widehat{\tilde{\mathbf{h}}}}_p^{(m,m)}\triangleq ({\mathbf{Z}}^{(m)}{)}^H{\mathbf{z}}_{:,p}^{(m)},{\mathbf{Z}}^{(m)}=\frac{({\mathbf{s}}^{(m)}\otimes {\mathbf{I}}_{D_m})}{\sqrt{E^{(m)}}{N}_{\text{F}}}$$

$${\widehat{\tilde{\mathbf{h}}}}_p^{(m,m)}={\tilde{\mathbf{h}}}_p^{(m,m)}+\sum _{\begin{array}{c}{m}^{\prime }=1\\ m^{\prime }\ne m\end{array}}^M\frac{\sqrt{E^{(m^{\prime })}}(({\mathbf{s}}^{(m)}{)}^H{\mathbf{s}}^{(m^{\prime })})}{\sqrt{E^{(m)}}{N}_{\text{F}}}{\tilde{\mathbf{h}}}_p^{(m,m^{\prime })}+{\mathbf{\xi }}_p^{(m)}$$

$${\mathbf{r}}_{:,p}=\left(\left[\sqrt{E^{(1)}}{\mathbf{s}}^{(1)}\dots \sqrt{E^{(M)}}{\mathbf{s}}^{(M)}\right]\otimes {\mathbf{I}}_R\right){\bar{\mathbf{h}}}_p^{(:)}+{\bar{\mathbf{\eta }}}_{:,p}$$

$${\widehat{\bar{\mathbf{h}}}}_p^{(:),\text{LS}}$$

$${\widehat{\bar{\mathbf{h}}}}_p^{(:),\text{MMSE}}$$

$${\mathbf{R}}_k^{(\text{ST},m)}\cong {\mathbf{R}}_p^{(m)},\text{for }p=(k-1)P+1,\dots ,kP.$$

$${\mathbf{x}}_{k+1}^{(\text{ST},m)}$$

$${\theta }_k^{(\text{ST},m)}$$

$${\widehat{\tilde{\mathbf{h}}}}_p^{(m,m)}\cong {\tilde{\mathbf{h}}}_p^{(m,m)}+{\mathbf{\xi }}_p^{(m)}$$

$$\mathbf{R}(\theta )\triangleq ({\mathbf{T}}^{(m)}{)}^H{\mathbf{R}}^{\text{FD}}(\theta ){\mathbf{T}}^{(m)},$$

$$\mathbf{R}(\theta )$$

$$\mathbf{E}(\theta )$$

$${\mathbf{f}}_k^{(m)}\triangleq {\left[\begin{array}{c}{\left[{\widehat{\tilde{\mathbf{h}}}}_{(k-1)P+1}^{(m,m)}\right]}^H\cdots {\left[{\widehat{\tilde{\mathbf{h}}}}_{kP}^{(m,m)}\right]}^H\end{array}\right]}^H.$$

$${\widehat{\tilde{\mathbf{h}}}}_p^{(m,m)}\cong \mathbf{E}({\theta }_k^{(\text{ST},m)}){\mathbf{\beta }}_p+{\mathbf{\xi }}_p^{(m)},$$

Full text

I Introduction

Massive multiple-input multiple-output (MIMO) is a widely studied technology that is the fundamental basis for new-generation communication systems [1]. Using large antenna arrays, massive MIMO increases the data rate, brings angular selectivity and angular focusing of the radiated power, and compensates for huge path loss in millimeter-wave communications [2]. It also enables spatial multiplexing of multiple users, removing the need for temporal or spectral division [3, 4]. Studies show that it increases the performance of full-duplex communications [5] and radar applications [6]. In addition, it is the basis for various new applications such as dual-functional radar-communication (DFRC) [7], reconfigurable intelligent surfaces (RIS) [8], and cell-free massive MIMO [9], and a myriad of 6G technologies such as high-speed train (HST) and vehicle-to-vehicle (V2V) communications [10].

As technology evolves, the number of actively communicating devices is expected to grow significantly [10]. In a time-varying scenario where multiple users are mobile, besides serving, channel estimation is a significant problem that results in explore/exploit dilemma and a huge computational load.

The hybrid beamformer (HBF) structure also complicates the estimation procedure. HBF is a widely offered structure for large antenna arrays instead of fully digital or fully analog beamformers for its efficiency regarding performance and physical viability [11]. It includes both analog beamformer (ABF) and digital beamformer (DBF) stages, with a reduced number of radio-frequency chains (RFCs). Since the ABF always needs to be directed to some angle(s), and angle of arrival (AoA) estimation is held in the digital stage after the ABF, AoA estimation is limited by the prior ABF beam coverage. In addition, user mobility might result in beam loss or unintended suppression due to ABF. These factors obligate either a scanning operation or a beam tracking procedure.

Computational complexity is a significant factor affecting the feasibility of a system design. In time-varying systems, varying channel parameters should be estimated frequently and systems using these parameters in their constructions should be updated accordingly. These estimation and update operations should be repeated within the interval in which the related parameters can be assumed to remain almost the same, namely the coherence time. On the other hand, the main task that brings the most complexity is interference mitigation, which requires inversion of matrices in the majority of techniques. Interference might be caused by the presence of multiple users and multipath components, or simply multiple signal clusters in a MIMO system.

For time-varying channels, the rate of variation imposed on channel parameters might be different [12]. Among the channel parameters, slowly varying and rapidly varying ones are often called slow-time (ST) and fast-time (FT) parameters, respectively. In the massive MIMO case, AoA and angular spread (AS) could be classified as ST parameters, while gains of individual micro-components inside a signal cluster could be FT parameters. Based on this classification, instantaneous channel estimation should be repeated after shorter intervals (in fast-time), while beam tracking can be repeated more rarely (in slow-time). Upon the estimates from these procedures, beamformer weights can be updated either in fast-time or slow-time depending on the method that is used, namely instantaneous or statistical beamforming.

I-A Literature Review

Angular estimation for massive MIMO systems is a widely studied topic. While exhaustive or hierarchical search is offered for initial estimation [13, 14, 15, 16, 17], lower-complexity beam tracking operations are offered after the initial phase to adaptively modify the estimate under variations through time.

In the vast majority of the studies, beam tracking is held after ABF consisting of phase shifters, which create discrete Fourier transform (DFT) beams, to probe angles in space. In [18, 19, 20, 21, 13], beam tracking means the selection of the optimum set of AoAs to be probed under time or physical (number of RFCs) limitations, after which the most powerful DFT beams are selected. The use of two perturbed beams around the main AoA is also common [22, 23, 24], from which the final AoA estimate is calculated accordingly. In the aforementioned studies, beam tracking is a phase where the beamformer is varied for a search. On the contrary, AoA variation is deduced from the variations on the current ABF output in [14, 25, 26] with the help of extended Kalman filter (EKF), and in [15] via maximum a posteriori (MAP) estimator. [27] and [15] employ fully digital beamformers, and [28] assumes erroneous AoA estimates and tracks the beams using recursive filtering.

Beam tracking in multiuser channels is investigated in only a few studies. Further, the near-far effect, to the authors’ knowledge, is not studied in the literature in the beam tracking context. That is the reason why the ABF output, whose sidelobe suppression is limited, is generally chosen for beam tracking. Against multiuser interference, [29] proposes a user scheduling algorithm, [27] employs orthogonal pilots, and [18] and [24] offer time or code domain multiplexing. In [20], multiuser estimation is performed either jointly by extending the MAP estimator, or sequentially by successive interference cancellation. [22] also applies joint estimation after ABF with the help of a maximum likelihood (ML) estimator.

Contrary to the algorithms that operate after ABF, [13] proposes measurement after a digital zero-forcing (ZF) beamformer after ABF, to suppress the multiuser interference further. However, since ZF is constructed with the channel estimates, it is, as we will define later, an instantaneous beamformer, which obligates frequent training phases. Indeed, beam tracking is generally mixed with the estimation of the rapidly changing channel gain coefficients in the literature. In [14, 25, 26, 22, 20], the proposed algorithms track the gain coefficient together with the AoA. On the contrary, [12] defines gain coherence time and angle coherence time such that the angle coherence time is much longer. Furthermore, [23] studies the optimum period for beam tracking.

Training sequence transmission is needed for beam tracking in most studies. Although uplink training is proposed generally, [19, 20, 21, 24] propose downlink training at the user equipment (UE), which requires feedback to the base station (BS). In [13], different DFT beams are tried inside the data mode, accepting data rate losses during the process.

Besides conventional communications, beam tracking is an important topic also for evolving technologies. It is studied for RIS [30], unmanned air vehicles (UAVs) [31], automated vehicles [32], and terahertz communications [33].

I-B Contributions

In this paper, a novel per-cluster estimation scheme is proposed for massive MIMO systems, which requires an inclusive consideration of the whole beamformer structure, beamforming approach, physical constraints, computational complexity, and time scheduling.

A statistical beamformer, namely generalized eigenbeamformer (GEB) from our previous work [34], is adapted to HBFs, which have a bank of analog phase shifters for ABF, to be used in the digital stage with the name reduced dimensional (RD) GEB (RD-GEB). Statistical beamformers do not depend on instantaneous channel coefficients, and they provide wider selection and null areas in the angular domain, based on the channel model that includes AS. Therefore, the validity of the beamformer lasts longer, which enables the definition of a separate coherent processing interval (CPI) for beam tracking, namely the slow-time CPI (ST-CPI), which is longer than the CPI for the instantaneous channel coefficients, fast-time CPI (FT-CPI). Moreover, it mitigates the multipath and multiuser interference better than widely proposed DFT beams, and creates interference-free subspaces for all clusters even under strong near-far effects. Consequently, multiuser channel estimation and beam tracking are performed via angular division in these subspaces in time-domain duplex (TDD) uplink mode, namely per-cluster, without spending resources via user scheduling algorithms, orthogonal pilot transmissions, or time division/synchronization between users. Further, the computational complexities of these estimation operations are very low because interference mitigation is handled previously by the statistical beamformer, and the dimension is further reduced. Also, beamformer construction and beam tracking are performed rarely in slow-time, but the instantaneous effective channel (IEC) estimator operates frequently with RD inputs in fast-time. Since the dominant complexity of the interference mitigation is undertaken by the ST beamformer, and the FT estimator is of low complexity, the time-averaged computational complexity also reduces. The indifference of the statistical beamformer to instantaneous variations of the channel is handled by a simple intra-cluster spatial channel matched filter (ICS-CMF), which is matched to the IEC.

The second major contribution is the design of novel beam tracking methods, namely beam-aware maximum likelihood estimator (BA-ML) and statistical extended Kalman filter (SEKF). Besides the fact that they both operate per-cluster after the statistical beamformer, both leverage the IEC estimates collected throughout an ST-CPI. Therefore, there is no need for an additional period of training. The difficulty of utilizing multiple channel estimates due to uncorrelated channel gains is handled meticulously in the design. While the BA-ML method is a variant of the nonlinear least squares method [35], SEKF is an EKF application whose observation is the second-order statistics, namely the effective channel covariance matrix (CCM). To sum up, the main contributions of this paper are listed below.

• •

The design of a novel per-cluster estimation approach with statistical beamforming, whose advantages are

• –

Suitability to

hybrid beamformers (HBFs),

• *

multiuser channels under notable near-far effect,

• *

time-varying channels and mobility,

• *

frequency-selective channels.

• –

High multiuser estimation performance thanks to the reduced dimensional interference-free subspaces, which also enables the usage of simpler estimators designed for single-user cases,

• –

Per-cluster estimation which removes the need for synchronization between the users,

• –

Low overall computational complexity thanks to the careful distribution of complexity load between FT and ST blocks, (especially the highly complex task of multiuser interference mitigation handled by ST statistical beamformer)

• –

Reduced delay spread for IEC estimation compared with joint estimation.

• •

The design of novel beam tracking methods, namely BA-ML and SEKF, which

• –

leverage the collection of IEC estimates,

• –

remove the need for an additional training period.

In addition, we adapt techniques from our prior work [34] to support the per-cluster estimation scheme for the remaining parts of the system, which are RD-GEB for beamformer construction, and beam-aware least squares (LS, BA-LS) for IEC estimation. Also, Orthogonal Matching Pursuit (OMP) method is adapted as a beam tracker to represent a joint estimation method in order to compare with per-cluster BA-ML and SEKF methods.

In the remainder of the paper, subscripts and superscripts generally indicate the time and cluster affiliations, respectively. Also, $\bar{x}$ and $\tilde{x}$ indicate the transform of an arbitrary variable $x$ after ABF and DBF, respectively. Also, $(\bm{X})_{a,b}$ indicates the entry of $\bm{X}$ in $a$th row and $b$th column. Finally, $\bm{x}^{H}$ and $\bm{x}^{T}$ indicate the Hermitian and transpose operations, respectively.

II System Model

This study considers a massive MIMO system where $U$ single-antenna mobile UEs are simultaneously communicating with a BS with a uniform linear array (ULA) of $N$ antenna elements in TDD uplink mode using single-carrier modulation. Signals arrive at BS as $M$ angularly resolved signal clusters, where $M~{}\geq~{}U$. Parameters of MIMO channels are grouped into two, namely FT and ST parameters. The BS estimates these parameters in separate repetition intervals, namely FT-CPI and ST-CPI. The estimated channel parameters can also be used in TDD downlink mode leveraging the channel reciprocity. The most important variables that will be defined in the subsequent sections are listed in Table I.

II-A Signal Model for Single Carrier Uplink Transmission

In TDD uplink mode, the received signal vector $\bm{y}_{n,p}\in\mathbb{C}^{N\times 1}$ in $n$th discrete time in $p$th FT-CPI is

[TABLE]

where $s_{n,p}^{(m)}$ are zero-mean unit-variance symbols, $E^{(m)}$ are cluster powers, and $\bm{h}_{p}^{(m)}\in\mathbb{C}^{N\times 1}$ is the channel vector for $m$th cluster. Also, $\bm{\eta}_{n,p}\sim\mathcal{CN}(\bm{0},N_{0}\bm{I}_{N})$ is the AWGN vector.

Signal clusters might belong to $U$ different users with $U=M$, or some clusters might be multipath components from a less number of users with $U<M$. This fact can be shown as

[TABLE]

where $b_{n,p}^{(u)}$ is the transmitted symbol from $u$th user, and $l_{m}$ indicates the delay for $m$th cluster. $\mathcal{U}(m)$ is the user-cluster association function from the set of clusters $\{1,\dots,M\}$ to the set of users $\{1,\dots,U\}$. It is simply an identity function $\mathcal{U}(m)=m$ for the case $U=M$.

II-B Multi-Cluster Massive MIMO Channel Model

In the $p$th FT-CPI, the channel model for the $m$th cluster is

[TABLE]

where $\bm{a}(\theta)\triangleq\frac{1}{\sqrt{N}}[1\,e^{j\pi\sin{\theta}}\cdots e^{j(N-1)\pi\sin{\theta}}]$, and $L$ is the number of rays. $\alpha_{l,p}^{(m)}\sim\mathcal{CN}(0,1)$ are complex gains which are i.i.d. in $l$ and $p$, and $\tilde{\theta}_{l,p}^{(m)}$ are azimuth angles of rays which are uniformly placed in the interval $(\theta_{p}^{(m)}-\Delta_{p}^{(m)}/2,\theta_{p}^{(m)}+\Delta_{p}^{(m)}/2)$. Therefore, a cluster has mean AoA $\theta_{p}^{(m)}$ and AS $\Delta_{p}^{(m)}$. As L goes to infinity, the channel model yields the cluster CCM $\bm{R}_{p}^{(m)}\triangleq\mathbb{E}\{\bm{h}_{p}^{(m)}(\bm{h}_{p}^{(m)})^{H}\}$ as

[TABLE]

at the $p$th FT-CPI. Therefore, from (1), the total covariance matrix $\bm{\Psi}_{p}\triangleq\mathbb{E}\{\bm{y}_{n,p}\bm{y}_{n,p}^{H}\}$ is expressed as

[TABLE]

II-C Time Variation Model

The block-fading model [36] is adopted to model the time variation in the channel. The channel $\bm{h}_{p}^{(m)}$ is assumed to be invariant within FT-CPIs of $T_{\text{F}}$ seconds, which consists of $N_{\text{F}}+N_{\text{S}}$ symbols as seen in Fig. 1. The channel varies after each FT-CPI as $p$ progresses.

In the channel definition given in (3), complex gains $\alpha_{l,p}^{(m)}$ are FT parameters, and mean AoA $\theta_{p}^{(m)}$ and AS $\Delta_{p}^{(m)}$ are ST parameters. Time variation is mainly governed by the uncorrelatedness of $\alpha_{l,p}^{(m)}$ in FT-CPI index $p$, while ST parameters also vary through $p$ in a slower and correlated manner.

In our work, variation of $\theta_{p}^{(m)}$ is modeled via a linear Gaussian state-space model while $\Delta_{p}^{(m)}$ is assumed to be constant. Transmitters are considered to move on a ring centered at the BS. $\theta_{p}^{(m)}$ and $\omega_{p}^{(m)}$ being the angular position and velocity in azimuth, respectively; the state vector in the $p$th FT-CPI is denoted by $\bm{x}_{p}^{(m)}\triangleq[\begin{matrix}\theta_{p}^{(m)}&\omega_{p}^{(m)}\end{matrix}]^{T}$. The state equation is

[TABLE]

where $\bm{\nu}_{p}^{(m)}\sim\mathcal{N}(\bm{0},\bm{\Sigma}_{\nu})$ is the innovation with $\bm{\Sigma}_{\nu}=\text{diag}\{\sigma_{\theta}^{2},\sigma_{\omega}^{2}\}$. The system initializes with a given $\bm{x}_{0}^{(m)}$.

II-D Practical Two-Stage Training Modes for Time-Varying Massive MIMO

The proposed structure of training and data transmission phases are illustrated in Fig. 1. FT-CPI and ST-CPI are defined to express the coherent interval in terms of FT and ST parameters. FT-CPIs consist of $N_{\text{F}}+N_{\text{S}}$ symbols, whereas ST-CPIs are longer and take $P$ FT-CPIs. In literature, they are also called gain and angle coherence time [12]. They determine the repetition time for the training modes FT training and ST beam tracking as described in Section IV, where beam tracking is performed more rarely. CPIs also determine the update times of different system blocks as described in Section III, where statistical blocks are updated more rarely.

This two-stage structure reduces the average computational complexity of beam tracking per time and enables low-complexity techniques for FT training with the aid of approximately known AoA, which is called beam-awareness. However, the design parameters should be chosen such that the duration for $N_{\text{S}}+N_{\text{F}}$ symbols and $P$ FT-CPIs are shorter than the actual FT-CPI and ST-CPI, respectively. As a design guideline, we provide these intervals in Table II for different hardware settings, where $f_{c}$, $W$, $v$, $d$, and $\lambda$ are carrier frequency, bandwidth, speed, distance, and wavelength, respectively.

FT-CPI, which corresponds to the channel coherence time [36], can be approximated as $0.1/f_{D}$, where $f_{D}=vf_{c}/c$ is the Doppler frequency considering the radial movement, and $c$ is the speed of light. Then, $0.1W/f_{D}$ gives the number of symbols in an FT-CPI. Therefore, there are approximately $c/(10vf_{c}/W)$ symbols in an FT-CPI. The time duration for an ST-CPI can be approximated as the beamwidth divided by angular speed. Beamwidth can be approximated as $2/N$ radians from $sin(\phi_{BW})\pi=2\pi/N$. Angular speed can be approximated as $v/d$ radians per second considering the tangential movement where $d$ is the distance between the BS and the receiver. Therefore, the number of symbols in an ST-CPI is $2dW/(Nv)$. Furthermore, the number of FT-CPIs in an ST-CPI, which is the limit for $P$, is $20(f_{c}d/c)/N$. Here, the ratio $f_{c}d/c=d/\lambda$ can be counted as a measure of required receive or transmit power, since it is directly related to the path-loss according to Friis transmission equation [36]. Therefore, a selection of $d/\lambda$ ratio, for example 10e3, represents a group of practices with similar power requirements, for example 100 meters and 30GHz, or 1 kilometer and 3 GHz.

III General System Structure

In this section, the general structure of the beamforming and processing system will be introduced, leaving the details of design procedures to Sections IV and V. The HBF structure is adopted as depicted in Fig. 2 for its practicality and energy efficiency. ABF produces inputs for $R$ RFCs via analog phase shifters and combiners, selecting all the interested clusters in the angular domain.

After the ABF, different than conventional systems, statistical beamforming is employed in cluster-specific DBF blocks via RD-GEB in Section V. DBFs create cluster-specific processing blocks, each of which includes an ICS-CMF and a per-cluster estimation block. The DBF is responsible for the formation of a cluster-specific subspace that mitigates the multi-cluster interference. Therefore, ICS-CMF and estimators are designed by ignoring the multiuser interference, yielding low complexity. ICS-CMF combines the outputs from the DBF as a matched filter. The chain of ABF, DBF, and ICS-CMF gradually reduces the signal dimension from $N$ to $R/M$, $D_{m}$, and 1, where $N>R/M\geq D_{m}\geq 1$. The output of DBF is also the place where the per-cluster estimation is performed, which consists of FT training and ST beam tracking, which will be detailed in Section IV.

Gains and phases of the outputs for clusters are corrected, and they are also combined by the digital cluster combiner (DCC) to reach symbol estimates for users, in the presence of multipath components ($U<M$) as described by (2). In this case, DCC is the block that brings the suitability for frequency-selective channels to the system.

In this structure, ABF and DBF are updated after each ST beam tracking rarely, while ICS-CMF and DCC are updated after each FT training more frequently.

III-A DFT-Based Analog Beamformer for Full Cluster Coverage

The analog beamforming matrix is denoted by $\bm{S}\in\mathbb{C}^{N\times R}$ for $R$ RFCs, whose columns are selected among DFT bases which results in a practical implementation via analog phase shifters. Also, the columns are orthonormal such that $\bm{S}^{H}\bm{S}=\bm{I}_{R}$. The output of the ABF is

[TABLE]

where $\bar{\bm{h}}_{p}^{(m)}\triangleq\bm{S}^{H}\bm{h}_{p}^{(m)}$ and $\bar{\bm{\eta}}_{n,p}\triangleq\bm{S}^{H}\bm{\eta}_{n,p}$. After the analog stage, covariance matrices $\bm{R}_{p}^{(m)}$ and $\bm{\Psi}_{p}$ are modified as

[TABLE]

where $\bar{\bm{\Psi}}_{p}\triangleq\mathbb{E}\{\bm{r}_{n,p}\bm{r}_{n,p}^{H}\}$ and $\bar{\bm{R}}_{p}^{(m)}\triangleq\mathbb{E}\{\bar{\bm{h}}_{p}^{(m)}(\bar{\bm{h}}_{p}^{(m)})^{H}\}$.

III-B Statistical Digital Beamformer for Inter-Cluster Separation

The digital stage starts with a bank of DBFs. DBFs are cluster-specific digital combiners for which statistical beamforming is employed. The DBF for the $m$th cluster $\bm{W}^{(m)}\in\mathbb{C}^{R\times D_{m}}$ has orthonormal columns such that $(\bm{W}^{(m)})^{H}\bm{W}^{(m)}=\bm{I}_{D_{m}}$. It processes the ABF output $\bm{r}_{n,p}$ and outputs $\bm{z}_{n,p}^{(m)}\triangleq\left(\bm{W}^{(m)}\right)^{H}\bm{r}_{n,p}$. Defining the total beamforming matrix $\bm{T}^{(m)}\triangleq\bm{S}\bm{W}^{(m)}$, the DBF output is

[TABLE]

where $\tilde{\bm{h}}_{p}^{(m,m^{\prime})}\triangleq(\bm{T}^{(m)})^{H}\bm{h}_{p}^{(m^{\prime})}$ is the IEC and $\tilde{\bm{\eta}}_{n,p}^{(m)}\triangleq(\bm{T}^{(m)})^{H}\bm{\eta}_{n,p}\sim\mathcal{CN}(\bm{0},N_{0}\bm{I}_{D_{m}})$. After the digital stage, covariance matrices of $\tilde{\bm{h}}_{p}^{(m,m^{\prime})}$ and $\bm{z}_{n,p}^{(m)}$, denoted by $\tilde{\bm{R}}_{p}^{(m,m^{\prime})}$ and $\tilde{\bm{\Psi}}_{p}^{(m)}$, respectively, are given in (11) and (12).

[TABLE]

III-C Intra-Cluster Spatial Channel Matched Filtering (ICS-CMF)

Different from temporal channel matched filtering, ICS-CMF compensates for the indifference of statistically constructed slow-time DBF to fast-time variations in the channel. After beamforming, $D_{m}$ outputs are obtained for the $m$th cluster. They are efficiently combined using the matched filtering method assuming the multi-cluster interference has already been suppressed to a negligible level by beamforming.111If multiple users ($\leq D_{m}$) fall into the angular sector of the same cluster, joint intra-cluster processing can be applied [34] both in data and training phases in RD subspace ($D_{m}$). The output of ICS-CMF is obtained by

[TABLE]

where $\hat{\tilde{\bm{h}}}_{p}^{(m,m)}$ is the estimate for the channel $\tilde{\bm{h}}_{p}^{(m,m)}$, which will be detailed in Section IV-A. The output is expressed as

[TABLE]

III-D Digital Cluster Combiner (DCC)

After ICS-CMF, the system has symbol estimates from different signal clusters as shown in (14). These symbol estimates have different gains and different delays. Furthermore, they might be multipath components from the same user, as described by (2). By substituting (2) into (14), we have

[TABLE]

To reach an estimate for the symbols $b_{n,p}^{(u)}$ from different users; delay, magnitude, and phases of symbol estimates from each cluster are aligned as a first step as

[TABLE]

which can be approximated as $\tilde{b}_{n,p}^{(m)}\cong b_{n,p}^{(\mathcal{U}(m))}+e_{n,p}^{(m)}$. Then, if the channel includes multipath components, implying $U<M$, symbol estimates for $M$ clusters are combined into symbol estimates for $U$ users, as

[TABLE]

for $u=1,\dots,U$, where $p_{m}^{(u)}$ is the combiner, which can be designed via various techniques, such as BLUE, maximal ratio combining, or selection combining. Note that the structure of $p_{m}^{(u)}$ is very simple such that $p_{m}^{(u)}\neq 0$ only for $m$ values for which $\mathcal{U}(m)=u$.

The digital cluster combiner is introduced to have an inclusive conceptual design. With this stage, the proposed scheme gains generality over frequency-selective and flat-fading channels. Also, it enables the per-cluster perspective of the proposed system design. After the advantages of the per-cluster scheme are exploited, per-user outputs are obtained via DCC in a simple way.

The design of $p_{m}^{(u)}$ is left out of scope, and its performance will not be evaluated since it requires a much more detailed scenario description. However, the quality of the symbol estimates at its input, which will be evaluated via SINR that will be defined in Section VII-B, determines its performance directly in any case.

IV Two-Stage Parameter Estimation Procedure

The task of channel estimation is held in two stages in the proposed system. They are named FT training and ST beam-tracking, which are introduced in Section II-D. To recite, there is an FT/ST separation in channel parameters, CPIs, estimation phases, and processing blocks in terms of their update rate.

The proposed per-cluster estimation approach employs separate estimators for each cluster as seen in Fig. 2. In Fig. 3, the estimation procedure is illustrated in more detail for one of the clusters. The fundamental aspect of the per-cluster estimation is that the estimation takes place after cluster-specific DBFs, which mitigate the interference better than conventional DFT beams. Also, the size of the effective channel decreases further. As a result, IEC estimators and beam trackers can be designed in a simpler way, as if only a single user exists, in the presence of multiple users or clusters.

On the other hand, in a conventional estimation scheme, the channel is estimated jointly after ABF, constructed with DFT beams. DFT beams are weak against interference, which charges the estimator with the interference mitigation task through pilots or other methods. Also, DBF depends on instantaneous channel estimates and needs more frequent updates, yet it is still responsible for interference mitigation for data transmission. Therefore, conventional methods suffer from performance loss, complexity increase, or a need for longer training overhead for channel estimation.

Beam trackers collect and use IEC estimates to estimate AoAs. It is very important that the channel estimator operates in fast-time and the beam tracker operates in slow-time, as shown by the dotted and dashed lines in Fig. 3. Therefore, ABF and DBF are updated in slow-time, while the ICS-CMF is updated in fast-time. Since DBF is slowly updated, a more complex design is tolerable which increases the interference mitigation performance.

IV-A FT Estimation: Instantaneous Effective Channel (IEC) Estimation

The IEC estimator aims to estimate the rapidly changing (instantaneous) RD channel observed after DBF, which is the reason why it is called the effective channel. Training sequences of length $N_{\text{F}}$ are transmitted, which are shown by the vector $\bm{s}^{(m)}=[s_{1+l_{m},p}^{(m)}\cdots s_{N_{\text{F}}+l_{m},p}^{(m)}]^{T}=[b_{1,p}^{(\mathcal{U}(m))}\cdots b_{N_{\text{F}},p}^{(\mathcal{U}(m))}]^{T}$, where $l_{m}$ is the cluster delay which is shown to be easy to track in Section IV-A3. The received signals from $n=1+l_{m}$ to $n=N_{\text{F}}+l_{m}$ are collected and the techniques below are applied. In the subsequent sections, $l_{m}=0$ is assumed for the sake of simplicity.

IV-A1 Beam-Aware Least Squares (BA-LS) Estimation

This technique observes the received signal after the ABF and DBF ($\bm{T}^{(m)}=\bm{S}\bm{W}^{(m)}$), which were set previously, to estimate the IEC denoted by $\tilde{\bm{h}}_{p}^{(m,m)}$. Beam-awareness implies that the subspace created by $\bm{T}^{(m)}$ is cluster-specific and interference mitigation is accomplished before by the beamformer, which yields a simple design. This technique also does not need synchronization between the sources since it operates per-cluster. Therefore, $\bm{s}^{(m)}$ for different $m$’s are uncorrelated.

The received signal $\bm{z}_{n,p}^{(m)}$ after $\bm{T}^{(m)}$ is as shown in (10). Received signals during the training are concatenated to obtain $\bm{z}_{:,p}^{(m)}\triangleq[(\bm{z}^{(m)}_{1,p})^{T}\dots(\bm{z}^{(m)}_{N_{\text{F}},p})^{T}]^{T}$, which is expressed as

[TABLE]

where $\tilde{\bm{\eta}}_{:,p}^{(m)}\sim\mathcal{CN}\left(\bm{0},\bm{I}_{N_{\text{F}}}\otimes(N_{0}\bm{I}_{D_{m}})\right)$, and $\otimes$ is the Kronecker product operator. Then, the received signal is processed by the LS estimator $\bm{Z}^{(m)}$ to obtain the LS estimate for the IEC $\tilde{\bm{h}}_{p}^{(m,m)}$ as

[TABLE]

where the expression for $\bm{Z}^{(m)}$ is simplified with the assumption of $(\bm{s}^{(m)})^{H}\bm{s}^{(m)}=N_{\text{F}}$. Consequently, the channel estimate can be expressed as

[TABLE]

where the second term is the inter-cluster interference and the last term is the error due to noise with $\bm{\xi}_{p}^{(m)}\sim\mathcal{CN}(\,\bm{0},\,\frac{N_{0}}{E^{(m)}N_{\text{F}}}\bm{I}_{D_{m}}\,)$.

IV-A2 Conventional Joint Estimation Techniques

In contrast to per-cluster BA-LS, a conventional estimation method might process the signals directly after the ABF and estimate the multi-cluster channels jointly. In this case, the $N_{\text{F}}R\times 1$ observation vector is $\bm{r}_{:,p}\triangleq\left[\bm{r}_{1,p}^{T}\cdots\bm{r}_{N_{\text{F}},p}^{T}\right]^{T}$ and the $MR\times 1$ regressor vector would be $\bar{\bm{h}}_{p}^{(:)}\triangleq[(\bar{\bm{h}}_{p}^{(1)})^{T}\cdots(\bar{\bm{h}}_{p}^{(M)})^{T}]^{T}$. The observation can be expressed as

[TABLE]

where $\bar{\bm{\eta}}_{:,p}$ is the AWGN with covariance $N_{0}\bm{I}_{N_{\text{F}}R}$. With these definitions and the ones in Section III-A, well-known LS and minimum mean squared error (MMSE) methods are implemented to estimate the concatenated channels $\bar{\bm{h}}_{p}^{(:)}$ jointly from the observation $\bm{r}_{:,p}$ after ABF [37]. The estimates are expressed as

[TABLE]

where $\bm{V}\in\mathbb{C}^{(N_{\text{F}}R\times MR)}$ is defined so that $\bm{r}_{:,p}=\bm{V}\bar{\bm{h}}_{p}^{(:)}+\bar{\bm{\eta}}_{:,p}$, and $\bar{\bm{R}}_{p}^{(:)}\in\mathbb{C}^{(MR\times MR)}$ is a block diagonal matrix which comprises $\bar{\bm{R}}_{p}^{(m)}$ for $m~{}=~{}1,\dots,M$. Note that these joint methods require synchronization between users before the training phase.

IV-A3 Discussion on Cluster Delay and Delay Spread

The discrete delay of $m$th cluster is denoted by $l_{m}$. The difference in delays of two clusters whose paths differ by $\Delta d$ is $\Delta l\triangleq|l_{m_{1}}-l_{m_{2}}|=\Delta d\times W/c$. For example, we have $\Delta l=20$ for $W=100$ MHz and $\Delta d=60$ m. Assuming constant normalized bandwidth $W/f_{c}$, $d\times W$ becomes directly related to path loss from Friis transmission equation [36]. Therefore, the range of a system with $W=10$ MHz increases compared to the one with $W=100$ MHz, and it could observe $\Delta d=600$ m, which results in $\Delta l=20$ again. This difference in delays actually increases the effective delay spread. The joint techniques should observe the received signal starting from the first symbol of the cluster with the minimum delay until the last symbol of the cluster with the maximum delay. Therefore, the observation vector in Section IV-A2 might be much longer in reality. On the contrary, the per-cluster estimator BA-LS observes only one of the clusters at a time, which removes this problem.

The estimation of the discrete delay is out of the scope of this paper. Once estimated, its tracking is straightforward. For example, for a source with speed $v$, consider its present position and that for $\Delta t$ later as two sources in the example above. Then, $\Delta d=v\Delta t$ in $\Delta l=\Delta d\times W/c$ expression. For $v=10$ m/s and $W=100$ MHz, $\Delta t$ should be 0.3 seconds to have $\Delta l=1$. That is $l_{m}$ changes by 1 after 0.3 seconds, which is a very long period of time compared with an FT-CPI, or even an ST-CPI. It can be tracked as a very-slow-time parameter.

IV-B ST Estimation: Beam Tracking

In the ST beam tracking phase, ST parameters of the channel, namely AoA $\theta_{p}^{(m)}$ and AS $\Delta_{p}^{(m)}$, should be estimated, which determine the second-order statistics of the channel as seen in (4). However, in this work, AS $\Delta_{p}^{(m)}$ is assumed to be constant and known.222AS estimation is studied in literature [12] and it is out of the scope of this paper. It is shown in [34] that the performance is not very sensitive to the exact value of AS. AoA of each cluster is estimated separately by cluster-specific estimators by leveraging the IEC estimates from the latest $P$ FT-CPIs. Therefore, beam tracking is actually a calculation phase, and transmission of a new training sequence is not needed. After the beam tracking phase, beamformers are updated with the new information.

IV-B1 Assumed Slow-Time Model

The overall design of ST operations, namely beam tracking and beamformer update, assume a coherent interval in terms of ST parameters, namely ST-CPI, which is longer than FT-CPIs as shown in Fig. 1. Therefore, these ST operations are repeated after a long time compared to FT operations. The length of ST-CPI is determined by the number of FT-CPIs within, denoted by $P$, as a design parameter, which should yield a similarity in terms of channel statistics such as $\bm{R}_{p}^{(m)}\cong\bm{R}_{p+P}^{(m)}$.

Due to the assumption of a different CPI, ST operations work under a mismatched channel model. Based on the aforementioned similarity, assumed variables related to previously defined ones are needed, indicated by the superscript ST. The first example of ST variables is

[TABLE]

As seen, ST variables are indexed by the ST-CPI index $k$, which is related to FT-CPI indices through $k=\lceil p/P\rceil$. Several ST variables will be used in the next sections without an explicit definition such as $\bm{\Psi}_{k}^{(\text{ST})}$, $\bar{\bm{R}}_{k}^{(\text{ST},m)}$, $\bar{\bm{\Psi}}_{k}^{(\text{ST})}$, $\bm{x}_{k}^{(\text{ST},m)}$, $\theta_{k}^{(\text{ST},m)}$, but they imply a similar relation.

The state-space model for angular variation defined in Section II-C is modified for ST operations as

[TABLE]

where $\bm{A}^{(\text{ST})}\triangleq\bm{A}^{P}$, $\bm{\nu}_{k}^{(\text{ST},m)}\sim\mathcal{N}(\bm{0},\bm{\Sigma}_{\nu}^{(\text{ST})})$, and $\bm{\Sigma}_{\nu}^{(\text{ST})}=\sum_{i=1}^{P}\bm{A}^{i}\bm{\Sigma}_{\nu}(\bm{A}^{i})^{T}$.

Furthermore, it is assumed that the overall beamformer structure has suppressed the multi-cluster interference to a negligible level. Therefore, the channel estimates expressed in (20) are suboptimally assumed as

[TABLE]

with $\hat{\tilde{\bm{h}}}_{p}^{(m,m)}\sim\mathcal{CN}(\bm{0},\bm{R}(\theta_{k}^{(\text{ST},m)})+\frac{N_{0}}{E^{(m)}N_{\text{F}}}\bm{I}_{D_{m}})$ in the design of estimators. The parametric RD CCM $\bm{R}(\theta)$ is expressed as

[TABLE]

where $\bm{R}^{\text{FD}}(\theta)\triangleq\frac{1}{\Delta}\int_{\theta-\Delta/2}^{\theta+\Delta/2}\bm{a}(\theta^{\prime})\bm{a}^{H}(\theta^{\prime})d\theta^{\prime}$ is the full dimensional (FD) CCM for mean AoA $\theta$ and AS $\Delta$. [28] shows that $\bm{R}^{\text{FD}}(\theta)\cong\,\text{diag}(\bm{a}(\theta))\bm{D}\,\text{diag}(\bm{a}(\theta))^{H}$ where $(\,\bm{D}\,)_{a,b}=\text{sinc}\left((a-b)\cos(\theta)\sin(\Delta/2)\right)$ for $a,b=1,\dots,N$. Further inspired by [28] for a simpler calculation of $\bm{R}(\theta)$, we can approximate it by

[TABLE]

where $\bm{E}^{\text{FD}}$ has $\sqrt{\lambda_{d}}\bm{e}_{d}$ in its columns where $\lambda_{d}$ and $\bm{e}_{d}$ are eigenvalues and eigenvectors of $\bm{D}$. Noting $\bm{R}^{\text{FD}}(\theta)$ and $\bm{D}$ are effectively low-rank for practical $\Delta$ values, $\bm{E}^{\text{FD}}$ can be constructed with the most dominant $D_{m}$ eigenvalues and eigenvectors, yielding a size of $N\times D_{m}$ where $D_{m}\ll N$.333The number of columns of $\bm{E}^{\text{FD}}$ and the number of DFT outputs are equal since both are determined by the effective rank of the CCM. $\bm{E}^{\text{FD}}$ is calculated and stored once, and $\bm{R}(\theta)$ can be calculated with $\mathcal{O}(ND_{m}^{2})$ multiplications via (29), instead of $\mathcal{O}(N^{2}D_{m})$ via (28), leveraging also the diagonal form of the $N\times N$ matrix $\text{diag}(\bm{a}(\theta))$ inside the $\bm{E}(\theta)$ expression.444$\bm{D}$ depends on $\theta$ through $\cos(\theta)\sin(\Delta/2)$, which is the transformed AS divided by $2\pi$, obtained from $\pi\sin(\theta+\Delta/2)-\pi\sin(\theta-\Delta/2)$ [28]. [34] shows the tolerance to AS errors, therefore dependency on $\theta$ can be neglected by selecting $\theta=0$. Nevertheless, one can store more than one $\bm{D}$ and $\bm{E}^{\text{FD}}$ matrices for quantized values of $\theta$ for a more accurate approximation.

IV-B2 Beam-Aware Maximum-Likelihood (BA-ML) Estimator

The observation vector used by the BA-ML estimator for the $m$th cluster in the $k$th ST beam tracking phase is

[TABLE]

In order to parameterize it on $\theta$, we rewrite (27) as

[TABLE]

where $\bm{\beta}_{p}\sim\mathcal{CN}(\bm{0},\bm{I}_{D_{m}})$ is the basis coefficient vector for the matrix $\bm{E}(\theta_{k}^{(\text{ST},m)})$ defined in (30), whose columns span the range space of $\bm{R}(\theta)$ in (28). The observation in (31) is parameterized on the sought unknown $\theta_{k}^{(\text{ST},m)}$ with the model in (32). However, $\bm{\beta}_{p}$ is also an unknown. In this case, the nonlinear least squares method, which is equivalently the ML method for Gaussian cases [35], yields the estimated $\theta_{k}^{(\text{ST},m)}$ from the observation $\bm{f}_{k}^{(m)}$ as

[TABLE]

which firstly finds and sets the maximizing set of $\{\bm{\beta}_{p}\}$ for each given $\theta$, then applies classical ML procedure on $\theta$. However, the given problem has $PD_{m}+1$ unknowns but $PD_{m}$ equations. Therefore, we reduce the size of $\bm{\beta}_{p}$ to $D_{m}^{\prime}<D_{m}$ and modify $\bm{E}(\theta_{k}^{(\text{ST},m)})\in\mathbb{C}^{D_{m}\times D_{m}}$ as $\bm{E}^{\prime}(\theta_{k}^{(\text{ST},m)})\in\mathbb{C}^{D_{m}\times D_{m}^{\prime}}$ by simply removing the weakest eigenvectors from $\bm{E}^{\text{FD}}$ in (30), which modifies (32) as

[TABLE]

Substituting (34) in (31), it is found in Appendix A that the AoA estimate can be calculated as

[TABLE]

where $\bm{M}(\theta)\triangleq\bm{I}_{D_{m}}-\bm{E}^{\prime}(\theta)\left(\bm{E}^{\prime H}(\theta)\bm{E}^{\prime}(\theta)\right)^{-1}\bm{E}^{\prime H}(\theta)$.

IV-B3 Statistical Extended Kalman Filter (SEKF)

A Bayesian method, particularly the Kalman filter, can exploit the correlation of AoA through time better due to slow variation. We propose the usage of second-order channel statistics as observation with EKF since it is a nonlinear function of AoA. The statistics are obtained via sample-mean covariance matrix, and the observation vector for EKF is obtained as

[TABLE]

where $\text{vec}\{\cdot\}$ is the vectorization operator. Assuming angular coherence, the observation can be rewritten as

[TABLE]

whose first term is the mean with $\bm{R}_{f}(\theta_{k}^{(\text{ST},m)})\triangleq\bm{R}(\theta_{k}^{(\text{ST},m)})+\frac{N_{0}}{E^{(m)}N_{\text{F}}}\bm{I}_{D_{m}}$, where (29) can be used for $\bm{R}(\theta)$. The second term $\bm{q}_{k}^{(m)}$ is the zero-mean error vector with covariance $\bm{Q}_{k}^{(m)}$, which is found in Appendix B as

[TABLE]

To implement EKF, $\text{vec}\{\bm{R}_{f}(\theta_{k}^{(\text{ST},m)})\}$ in (37) is linearized using the first-order Taylor expansion, which yields

[TABLE]

where $\bm{B}_{k}^{(m)}$ is the Jacobian matrix of $\text{vec}\{\bm{R}_{f}(\theta)\}$ with respect to $\bm{x}_{k}^{(\text{ST},m)}$, to be calculated at $\hat{\bm{x}}_{k|k-1}^{(\text{ST},m)}=[\hat{\theta}_{k|k-1}^{(\text{ST},m)}\quad\hat{\omega}_{k|k-1}^{(\text{ST},m)}]$. It is found as

[TABLE]

whose second column, derivative with respect to speed, is zero due to the angular coherence assumption. The derivative with respect to the angular position can be calculated numerically or analytically using the simple structure in (29).

Using these definitions, well-known iterations of EKF [37] are applied for each cluster $m$. These iterations are shown below for the $k$th beam tracking phase.

[TABLE]

[TABLE]

After the $k$th ST beam tracking, beamformers are updated with predicted AoAs $\hat{\theta}_{k+1|k}^{(\text{ST},m)}$ for $m=1,\dots,M$, which are the first elements of $\hat{\bm{x}}_{k+1|k}^{(\text{ST},m)}$ for $m=1,\dots,M$.

IV-B4 Orthogonal Matching Pursuit (OMP)

Previous methods work per-cluster in the subspace created by the DBF $\bm{W}^{(m)}$, under our proposed per-cluster estimation scheme. Alternatively, conventional methods might directly work on the output of the ABF $\bm{S}$. We adapt the well-known OMP technique to our framework to obtain this alternative.

OMP is a compressed sensing method that solves the systems in the form of $\bm{f}=\bm{G}\bm{x}$, where $\bm{G}$ is a fat matrix. For angular estimation, $\bm{f}$ and $\bm{G}$ could be set as the received signal vector and a matrix with steering vectors at the angles to be searched in its columns, respectively. Then, the angular estimates could be the angles of the columns of $\bm{G}$ related to the dominant values in $\bm{x}$. However, this form could fail when the observation is noisy and weights are zero-mean random variables, that is, the channel is subject to fading.

The aforementioned adaptation is about the training sequence usage, presence of a preprocessing (ABF), countermeasures to the near-far effect, and most importantly, the utilization of the previous $P$ estimates. The challenge in the last factor is that the channel gains are uncorrelated in the previous $P$ FT-CPIs, and careless integration of them might result in the weights averaging out to zero.

For our adaptation, the system to be solved is $\bm{F}=\bm{G}\bm{X}$ where $\bm{F}\in\mathbb{C}^{RN_{\text{F}}\times P}$ includes all the observed signals at the FT-training phases from the previous $P$ FT-CPIs. For the $k$th beam tracking phase, it can be written as

[TABLE]

where $\bm{r}_{:,p}$ is as given in Section IV-A2. The matrix $\bm{G}\in\mathbb{C}^{RN_{\text{F}}\times MN_{\theta}}$ is in the form of $\bm{G}\triangleq\left[\bm{G}^{(1)}\cdots\bm{G}^{(M)}\right]$ with

[TABLE]

where $N_{\theta}$ is the number of angles to be searched per cluster. The $c$th column of $\bm{G}$, denoted by $\bm{G}_{:c}$, is affiliated with a cluster and an AoA. Let the functions $M^{OMP}(c)$ and $\theta^{OMP}(c)$ map these columns to clusters and AoAs, respectively. With these definitions, the modified OMP algorithm is given in Algorithm 1.

Considering the system $\bm{F}=\bm{G}\bm{X}$ with these definitions, it is seen that the rows of the weight matrix $\bm{X}\in\mathbb{C}^{MN_{\theta}\times P}$ are actually an angular map of signal presence monitored throughout $P$ FT-CPIs. In this regard, it is expected to be row-sparse, that is, only a few rows have significant nonzero content. Therefore, we adapt the selection mechanism of the columns of the matrix $\bm{G}$ in the conventional OMP method as in step 6 of the Algorithm to reflect the row-sparsity of the matrix $\bm{X}$. The proposed mechanism also solves the problem of coherency, where the channel gains are uncorrelated through FT-CPIs. In step 7, estimates are categorized according to cluster affiliation, which is a general problem in joint estimation techniques. The categorization is, in fact, strengthened by the difference of the training sequences as seen in (47). Then, steps from 8 to 11 handle the multiple estimates from the same cluster, which is needed due to AS and near-far effect. Note that this method requires pilot transmission and therefore synchronization between users, unlike the proposed per-cluster methods BA-ML and SEKF.

V Statistical Beamformer Construction

V-A Full Dimensional GEB (FD-GEB)

GEB is a statistical beamformer whose performance is widely evaluated in [34]. Similar to the well-known Capon beamformer [35], GEB minimizes interference while the intended signal power is kept constant. However, GEB accomplishes this task over variances as

[TABLE]

where $\bm{w}$ is a beamformer vector, and the terms are variances of $\bm{w}^{H}\bm{y}_{n,p}$ and $\bm{w}^{H}\bm{h}_{p}^{(m)}$, respectively. The covariance matrices can be calculated parametrically, and $c$ is a constant. The solution to this problem requires generalized eigendecomposition of the matrix pair $(\bm{R}_{k}^{(\text{ST},m)},\bm{\Psi}_{k}^{(\text{ST})})$, expressed as

[TABLE]

Choosing $\bm{w}_{\text{opt}}=\bm{e}_{\text{max}}$ solves the problem in (48) and yields the maximum signal-to-interference-plus-noise ratio (SINR), where the generalized eigenvector $\bm{e}_{\text{max}}$ corresponds to the maximum generalized eigenvalue $\lambda_{max}$.

FD-GEB is a single-stage beamformer where $\bm{T}^{(m)}\in\mathbb{C}^{N\times D_{m}}$ is constructed with the most dominant $D_{m}$ generalized eigenvectors, using either a fully digital or a fully analog beamformer (with the double phase shifter structure).

V-B Reduced Dimensional GEB (RD-GEB)

For HBFs, GEB can be implemented in the subspace created by an ABF $\bm{S}$ using the effective CCM pair $(\bar{\bm{R}}_{k}^{(\text{ST},m)},\bar{\bm{\Psi}}_{k}^{(\text{ST})})$, given in (8) and (9), instead of full-dimensional CCM pair $(\bm{R}_{k}^{(\text{ST},m)},\bm{\Psi}_{k}^{(\text{ST})})$. However, we propose a slightly different technique that is more robust and integrated with the design of the ABF $\bm{S}$. After the AoA estimates $\hat{\theta}_{k}^{(\text{ST},m)}$ for $m=1,\dots,M$ are produced in $k-1$st ST beam tracking phase, ABF is formed by the selection of the $R$ DFT bases which are directed to the closest angles to $\hat{\theta}_{k}^{(\text{ST},m)}$. The DFT frequencies are

[TABLE]

and the $k$th DFT basis is $\bm{u}(\phi_{k})$ where $\bm{u}(\phi)\triangleq\frac{1}{\sqrt{N}}[\begin{matrix}1&e^{j\phi}&\dots&e^{j(N-1)\phi}\end{matrix}]^{T}$. Then, $R/M$ DFT frequencies are selected for each $m$ such that $\sum_{r=1}^{R/M}|\phi_{k_{r}^{(m)}}-\pi\sin(\hat{\theta}_{k}^{(\text{ST},m)})|$ is minimized, where $k_{r}^{(m)}$ are the indices for the selected set.555Note that $k_{r}^{(m)}$ are consecutive in $r$ for an arbitrary $m$, that is $k_{r}^{(m)}~{}=~{}\bar{k}^{(m)}+r$ for some $\bar{k}^{(m)}$. Therefore, the selection is straightforward. Then the ABF $\bm{S}$ is constructed as

[TABLE]

We assume that each possible DFT basis selects an angular sector, such that $\bm{u}(\phi_{k})$ selects the angular sector from $\theta_{1}$ to $\theta_{2}$, where $\pi\sin(\theta_{1})=\phi_{k}-\pi/N$ and $\pi\sin(\theta_{2})=\phi_{k}+\pi/N$. These angular sectors form a codebook of covariance matrices $\{\bm{C}_{k}\}$ for $k=1,\dots,N$, where

[TABLE]

Using this codebook, approximated CCMs for the clusters are calculated as

[TABLE]

Then, the DBF $\bm{W}^{(m)}$ is obtained using the most dominant $D_{m}$ generalized eigenvectors of the matrix pair $(\hat{\bar{\bm{R}}}_{k}^{(\text{ST},m)},\hat{\bar{\bm{\Psi}}}_{k}^{(\text{ST})})$, which yields the total beamformer as $\bm{T}^{(m)}=\bm{S}\bm{W}^{(m)}$.

Columns of $\bm{W}^{(m)}$ (generalized eigenvectors) are not necessarily orthogonal, therefore they are orthogonalized via QR decomposition, assigning the obtained vectors as the new $\bm{W}^{(m)}$ such that $(\bm{W}^{(m)})^{H}\bm{W}^{(m)}=(\bm{T}^{(m)})^{H}\bm{T}^{(m)}=\bm{I}_{D_{m}}$. This process does not bring any information loss and the range space of $\bm{W}^{(m)}$ is preserved. In addition, it has less computational complexity than eigendecomposition, and it provides numerical stability and notational simplicity.

VI Computational Complexity

Computational complexities of all the alternatives are tabulated in Table III in terms of the number of multiplications using $\mathcal{O}(\cdot)$ notation. The dominant operations resulting in these complexities are also given in the third column where MI, MM, and ED mean matrix inversion, matrix multiplication, and eigendecomposition, respectively. Also * means a size reduction in inversion is possible and applied in favor of related techniques. Proposed per-cluster operations BA-LS, BA-ML, SEKF, and GEB construction are additionally multiplied by $M$ since they are applied for each cluster separately. In addition, all beam tracking methods and GEB construction are applied once in every $P$ FT-CPIs. Instantaneous channel estimators include the training size $N_{\text{F}}$, and search-based beam trackers BA-ML and OMP include $N_{\theta}$, which is the number of AoAs to be searched per cluster. For joint LS and MMSE channel estimators, and the MMSE beamformer, reduced matrix inversion sizes are shown using properties of Kronecker product and Woodbury matrix identity, although they are originally larger. Multiplications with matrices expressed by a Kronecker product, especially when it includes an identity matrix, might also yield a reduced number of multiplications when inspected in detail, such as the multiplication in (19), which yields $\mathcal{O}(MN_{\text{F}}D_{m})$ instead of $\mathcal{O}(MN_{\text{F}}D_{m}^{2})$.

Joint LS channel estimation in (22) has the inverse of a matrix of size $MR\times MR$ and multiplication of matrices of sizes $MR\times MR$, $MR\times N_{\text{F}}R$, and $N_{\text{F}}R\times 1$. Joint MMSE channel estimation in (23) has the inverse of a matrix of size $N_{\text{F}}R\times N_{\text{F}}R$ and multiplication of matrices of sizes $MR\times N_{\text{F}}R$, $N_{\text{F}}R\times N_{\text{F}}R$, and $N_{\text{F}}R\times 1$. The size of the matrix inversion can be reduced to $M$ for Joint LS using a property of Kronecker product, and to $MR$ for Joint MMSE using the Woodbury matrix identity. On the other hand, MMSE BF construction in (56) is dominated by a matrix inversion of size $R\times R$ or $N\times N$, however, it can be reduced to an $M\times M$ inverse by the Woodbury matrix identity. Usage of these properties and identities changes also the placement of matrices and matrix multiplication strategy, and the reduced complexities in Table III are obtained. For GEB, eigendecomposition for $x\times x$ matrices has a complexity of $\mathcal{O}(x^{3})$, where this size is $N$ for FD-GEB and $R$ for RD-GEB.

We can roughly assume that $M\cong D_{m}\cong R/M<N_{\text{F}}\ll N_{\theta}<N<P$ (see Table IV). Accordingly, we can sort the computational complexity of techniques as BA-LS < Joint LS $\ll$ Joint MMSE for instantaneous channel estimation, SEKF < BA-ML $\ll$ OMP for beam tracking, and RD-GEB $\ll$ FD/RD MMSE BF < FD-GEB for beamformer construction.

Consequently, it is clearly seen that the proposed per-cluster estimation approach with statistical beamforming, represented by BA-LS, BA-ML, SEKF, and GEB, yields very low computational complexity compared with the conventional joint estimation approach with instantaneous beamforming, represented by joint LS, joint MMSE, OMP, and MMSE BF.

VII Numerical Evaluations

VII-A Auxiliary Definitions

For comparison purposes in numerical results, some auxiliary definitions are needed. Firstly, the prefix RD (reduced dimensional) implies that the DBF is working after an ABF in an HBF system, whereas FD (full dimensional, or fully digital) refers to a system where the beamformer is fully digital and there is no ABF, therefore $\bm{S}=\bm{I}_{N}$ and there are $R=N$ RFCs.

DFT beamformer has the same ABF $\bm{S}$ as RD-GEB. However, the DBF $\bm{W}^{(m)}$ is a matrix of ones and zeros, only selecting the RFCs (columns of $\bm{S}$) for the $m$th cluster in (51) so that $\bm{T}^{(m)}=\bm{S}\bm{W}^{(m)}=\tilde{\bm{S}}^{(m)}$. This setting of DBF can be seen as the most primitive attempt to create cluster-specific subspaces after ABF. The improvement brought by RD-GEB can be seen from its difference from DFT BF.

The periodogram is a spectral estimation method [35], which is the implicit basis for a variety of techniques in the beam tracking literature. It is applied in the range space of the DFT beamformer in the case of HBF structures. The peak of $\rho(\theta)\triangleq\sum_{n,p}|\bm{a}(\theta)^{H}\tilde{\bm{S}}^{(m)}\bm{z}_{n,p}^{(m)}|^{2}$ is searched for the $m$th cluster, where $\bm{z}_{n,p}^{(m)}=(\tilde{\bm{S}}^{(m)})^{H}\bm{y}_{n,p}$. Inspecting (52), note that $\rho(\theta)$ measures the power at the RFC outputs one by one at certain $\theta$ values. Therefore, it represents the related studies in Section I-A. Considering all the training phases in an ST-CPI are used, the summation above is taken from $PN_{\text{F}}$ samples.

MMSE beamformer is another beamforming method, different than MMSE channel estimator, for which an MMSE estimator [37] is constructed where symbols $s_{n,p}^{(m)}$ for $m=1,\dots,M$ are treated as the multiple unknown parameters with the observation $\bm{y}_{n,p}$ in (1) or $\bm{r}_{n,p}$ in (7), depending on whether it is FD or RD. For the RD case,

[TABLE]

where $\bar{\bm{H}}_{p}\triangleq[\bar{\bm{h}}_{p}^{(1)}\dots\bar{\bm{h}}_{p}^{(M)}]$. As seen, the channels $\bm{h}_{p}^{(m)}$ or $\bar{\bm{h}}_{p}^{(m)}$ are needed for beamformer construction because they are combiners of unknown parameters and should be known. That is why we categorize the MMSE beamformer as an instantaneous beamformer since it should be updated with each new instantaneous channel estimate in the rate of channel decorrelation, that is FT-CPI, after each FT training.

There are two modes related to beam tracking, the Self-Driven mode and the Genie-Aided BF mode. In the proposed system, IEC estimators and beam trackers operate after the beamformers, which are constructed via the previous AoA estimates. Therefore, all the performance measures depend on the previous AoA estimate which describes the Self-Driven mode. On the other hand, the Genie-Aided BF mode removes the effect of the previous AoA estimate for performance analysis purposes. The beamformers are constructed with true AoAs at the beginning of each ST-CPI, and all the performance measures are collected at the end of each ST-CPI, including angular RMSE from beam trackers although their AoA estimates will not be used in the next ST-CPI. Note that this does not mean that the mismatches due to the movement inside the ST-CPIs are eliminated.

VII-B Performance Measures

VII-B1 Average NMSE for BA-LS

The performance of IEC estimation is measured by normalized mean squared error (NMSE), which is defined for the $m$th cluster as666Note that the expectations both in SINR and NMSE calculations are FT expectations that are conditioned on the settings at the related FT-CPI, such as beamformers and positions.

[TABLE]

It is shown in Appendix C that the NMSE can be calculated as

[TABLE]

Then, the average NMSE is found by averaging $\text{NMSE}_{p}^{(m)}$ through all the possible $p$ instants in Monte-Carlo experiments.

VII-B2 Angular RMSE for Beam Tracking

Angular error is defined as $\hat{\theta}_{k}^{(\text{ST},m)}-\theta_{(k-1)P+1}^{(m)}$, and the angular root mean squared error (RMSE) is

[TABLE]

for the $m$th cluster, where $K$ is the number of all the possible slow-time instants in Monte-Carlo experiments. Note that the error definition neglects the variation inside the ST-CPI and focuses only on the estimation performance.

VII-B3 Average SINR after ICS-CMF

The overall performance of the ABF, DBF, beam tracker, and IEC estimator will be measured via SINR at the output of ICS-CMF, which is the symbol estimate $\hat{s}_{n,p}^{(m)}$ given in (14). Since the estimated channel is known to the detector rather than the true channel, the true signal term is taken as $S_{p}^{(m)}\triangleq\sqrt{E^{(m)}}(\hat{\tilde{\bm{h}}}_{p}^{(m,m)})^{H}\hat{\tilde{\bm{h}}}_{p}^{(m,m)}s_{n,p}^{(m)}$, and the interference-plus-noise term $N_{p}^{(m)}=\hat{s}_{n,p}^{(m)}-S^{(m)}$ is the remaining part. With these definitions, SINR for the $m$th cluster is defined as

[TABLE]

Then, the average SINR is found by averaging $\text{SINR}_{p}^{(m)}$ through all the possible $p$ instants in Monte-Carlo experiments.

VII-C Simulation Settings

Selected parameters for simulations are listed in Table IV. One FT-CPI consists of $N_{\text{F}}+N_{\text{S}}=1000$ symbols, which is suitable according to Table II, and it takes $T_{\text{F}}=10\mu s$ assuming a bandwidth of 100 MHz. Selection of $N_{\text{F}}=10$ and $N_{\text{S}}=990$ yields 1% training overhead.

There are $M=4$ signal clusters. The first cluster, which will be the most interested one, suffers from a significant near-far effect since it is surrounded by 20 to 30 dB stronger clusters. The movement model given in Section II-C is implemented with $\sigma_{\theta}^{2}=1.45e-4$, $\sigma_{\omega}^{2}=1.46e-6$, and initial AoAs given in Table IV. This setting of variances corresponds to 10 m/s speed and 1 m/s2 acceleration assuming 150 meters distance from the BS after 1 second of random movement.777We interpret the standard deviation for the angular position (speed) divided by time as speed (acceleration), in the innovation covariance matrix after 100.000 FT-CPIs (1 second), which is computed similarly to $\bm{\Sigma}_{\nu}^{(\text{ST})}$ in Section IV-B1. One example of the resultant movement process is plotted in Fig. 4 for a duration of 1 second.888Note that the selection of $T_{\text{F}}=10\mu s$ with the assumption of 100 MHz bandwidth is just to give an insight on a practical operation. $T_{\text{F}}$ only affects the movement model, and variances would be changed for different $T_{\text{F}}$ and bandwidth just to have a similar movement with the one in Fig. 4.

Parameters $D_{m}$ and $R/M$ should be comparable with the practical rank of effective CCMs $\tilde{\bm{R}}_{p}^{(m,m)}$, which is determined by the AS $\Delta^{(m)}$ and the number of antennas $N$. If the ratio $R/N$ is limited, the number of active users can be reduced by user scheduling, or $D_{m}$ and $R/M$ can be reduced by treating the AS $\Delta^{(m)}$ as a design parameter.

For each cluster, only the selected region by ABF $\tilde{\bm{S}}^{(m)}$, which is discussed below (51), is searched for BA-ML, OMP, and the periodogram with a resolution of $0.1^{\circ}$, which yields $N_{\theta}\cong 40$. $D_{m}^{\prime}=2$ is chosen for SEKF.

VII-D Numerical Results

In this section, numerical results from the simulations will be shared. In the simulations, $p_{\text{max}}$ FT-CPIs of movement period, whose one realization for $p_{\text{max}}=1e5$ is given in Fig. 4, is repeated many times to have a Monte-Carlo experiment.999If any two clusters get closer than $3^{\circ}$, or any cluster gets outside of [$-60^{\circ}$,$60^{\circ}$] region, related realization is ended before $p_{\text{max}}$ FT-CPIs. In the simulations, the parameters given in Section VII-C are used, and the performance measures given in Section VII-B are collected.

As a beginning, the motivation behind the proposed per-cluster estimation approach with statistical beamforming will be clarified via comparisons with conventional techniques in Figures 5 and 6. In these figures, mobility and beam tracking is temporarily left out of focus with the settings of Genie-Aided BF mode and $p_{\text{max}}=1000$, for which mobility is almost absent as seen in Fig. 4. In Fig. 5, statistical beamforming (FD-GEB, RD-GEB, and DFT BF) and instantaneous beamforming (FD/RD MMSE BF) are compared according to beamformer update time $P$, which work with second-order statistics (CCMs) and actual channel estimates, respectively. As expected, slow-time alternatives are robust against $P$, while MMSE beamformers fail immediately when $P>1$, since the actual channel decorrelates but statistics remain almost the same. Although the best performance is reached by MMSE beamformers using MMSE channel estimates with $P=1$, the optimum choice in terms of computational complexity per time is statistical beamformers with a large $P$. The first scheme needs MMSE channel estimates for each FT-CPI which is computationally complex as seen from Table III, and the second scheme needs channel estimates from BA-LS, which is simpler. In addition, RD-GEB has less averaged complexity than the RD MMSE beamformer. The simpler alternative of channel estimation for instantaneous beamforming, joint LS, causes a performance loss. Finally, the RD alternatives, which are suitable for HBFs, exhibit minimal performance loss compared with FD alternatives, which are structurally and computationally more complex and needy.

In Fig. 6, instantaneous channel estimation methods are compared against signal-to-noise ratio (SNR). As discussed in Section VI, BA-LS is computationally simpler than the conventional techniques of joint LS and joint MMSE. Although simpler, it is seen in Fig. 6 that BA-LS does not lose performance compared with the superior technique joint MMSE.101010The vertical axis is normalized MSE, and the reason for BA-LS to seem unexpectedly better than joint MMSE is the difference in the sizes of the channels, which are $D_{m},D_{m},R/M,N,N,R,R$ in the order of the legend. This high performance of BA-LS arises from the prior processing via DBF, where GEB is employed. The role of GEB is seen from the performance loss of BA-LS with DFT BF, compared to BA-LS with RD-GEB case. In conclusion, BA-LS with GEB is a very efficient technique in terms of performance and computational complexity.

The self-driven mode, which is defined in Section VII-A might suffer from instability. Fig. 7 depicts the angular RMSE when the beamformers are constructed with a given error in the AoA. In this regard, Fig. 7 constitutes a basis for the stability of self-driven operation. For a stable self-driven performance, the next angular error is desired to be smaller than the given error, which causes a given error to gradually get smaller and become zero. The region where this behavior is expected is called the lock-in region, as shown in the figure. First of all, it is seen that the error in the next estimate is step-like against the previously given error, which is due to discrete angles in the definition of ABF and discretized angular sectors in the definition of RD-GEB. Second, note that the AS is 3 degrees, and $R/M=4$ RFCs of ABF create a beamwidth of approximately 3.5 degrees. Therefore a given error larger than 3 or less than -3 means total beam loss. It is seen in Fig. 7 that BA-ML and OMP offer the most robust performance when there is no beam loss. SEKF also remains in the lock-in region, however, the error due to linearization in (39) and (40) increases with the given angular error. Periodogram fails due to the near-far effect because the cluster $m=2$ is at 20 degrees, which is 30 dB stronger, and estimates deviate towards this cluster. Lastly, we see that the OMP shows nearly 0.5 degrees of RMSE when the given angular error is zero. This is because OMP actually searches for a point target, but BA-ML and SEKF search for a cluster with some AS through their modeling by (34) and (37). Therefore, OMP appoints the first powerful AoA candidate inside the AS as an estimate, but others find the center of the observed cluster.

In Figures 8, 9, and 10, performances of various angular estimators, the proposed IEC estimator BA-LS, and overall performance measure SINR will be shared. Although some results for the Genie-Aided BF mode are also given for comparison, the main focus is on the performance of the Self-Driven mode with the setting of $p_{\text{max}}=100000$, which leads to more angular variation. In Fig. 8, tracking behaviors of different angular estimators are given. This figure shows more clearly that the actual AoAs are changing after each FT-CPI, but estimates are produced after every $P$ FT-CPIs, or an ST-CPI. Inside an ST-CPI, beamformers $\bm{S}$ and $\bm{W}^{(m)}$, which are constructed upon these estimates, remain unchanged. It is also seen that the periodogram immediately deviates towards the stronger neighboring cluster, while others are robust to the near-far effect.

In Fig. 9, cumulative distribution functions (CDFs) of angular error for different beam trackers are given. It is seen that the proposed estimators SEKF and BA-ML outperform the alternatives with 90% of estimates being smaller than 0.5 degrees when operated with FD or RD-GEB. Comparing these with OMP after ABF, we see the contribution of applying cluster-specific DBFs before the estimation, that is the per-cluster approach. With this approach, the complexity is decreased and performance is enhanced.

In Fig. 10, the CDF of NMSE of the proposed instantaneous channel estimation technique BA-LS is plotted. In this figure, the effects of various beamforming and beam tracking techniques on the instantaneous channel estimation accuracy are shown. It is seen that the GEB alternatives used together with the proposed BA-ML and SEKF beam trackers, which are the representatives of the per-cluster estimation approach, result in NMSE values between 0.035 and 0.05. Note that these numbers are lower-bounded by $(\text{SNR}\times N_{\text{F}})^{-1}$, which is 0.01 for $m=1$ and $N_{\text{F}}=10$. It is seen that RD-GEB works with minimal performance loss compared with FD-GEB, and proposed BA-ML and SEKF beam trackers yield NMSE performance similar to Genie-Aided BF mode. Finally, the huge advantage of per-cluster operation is seen from the difference between RD-GEB and DFT BF, which are different only in the design of $\bm{W}^{(m)}$.

As a final performance measure, Fig. 11 shows angular RMSE and SINR performances against the design parameter $P$, the number of FT-CPIs in an ST-CPI. There is a performance loss at large $P$ for all the alternatives due to the loss of angular coherency because of mobility. For small $P$, BA-ML performance degrades due to the decrease in the observation period. However, it is seen that SEKF is not affected by this factor since it is a Kalman filter variant that can collect information from the previous ST-CPIs. Therefore, SEKF can succeed even with the primitive DFT BF, but DFT BF fails in beamforming and channel estimation as seen in subplot (b) and also in Fig. 10. Finally, it is seen that the BA-ML and SEKF methods in Self-Driven mode perform similarly to the Genie-Aided mode. Also, the RD alternatives with hybrid beamforming perform with minimal loss compared with the structurally complex and needy FD alternatives.

VIII Conclusion

In this paper, a novel channel estimation framework is proposed, which includes IEC estimation and beam tracking. This framework, namely the per-cluster estimation with statistical beamforming, is based on the spatial multiplexing of users and signal clusters both in data and training periods. Also, instantaneous channel estimation is repeated in fast-time while beam tracking and beamformer update are slow-time operations. The proposed scheme is shown to be advantageous in terms of both computational complexity and performance compared with the conventional approaches.

Appendix A Angular Estimate for BA-ML

In this appendix, (35) will be derived using (33) and the definitions in IV-B2. Since the channel estimates in (31) are taken from different FT-CPIs, they are independent and we can use $p(\bm{f}_{k}^{(m)}|\theta,\{\bm{\beta}_{p}^{\prime}\})=\prod_{p=1}^{P}p(\hat{\tilde{\bm{h}}}_{p}^{(m,m)}|\theta,\bm{\beta}_{p}^{\prime})$. Then,

[TABLE]

We assume that $p(\hat{\tilde{\bm{h}}}_{p}^{(m,m)}|\theta,\bm{\beta}_{p}^{\prime})$ is expressed by $\mathcal{CN}(\bm{E}^{\prime}(\theta)\bm{\beta}_{p}^{\prime},\frac{N_{0}}{E^{(m)}N_{\text{F}}}\bm{I}_{D_{m}})$. Then, the likelihood $p(\hat{\tilde{\bm{h}}}_{p}^{(m,m)}|\theta,\bm{\beta}_{p}^{\prime})$ is maximized for given $\hat{\tilde{\bm{h}}}_{p}^{(m,m)}$ and $\theta$ by the LS estimate $\hat{\bm{\beta}}_{p}$ of $\bm{\beta}_{p}^{\prime}$. Then, the mean $\bm{E}^{\prime}(\theta)\hat{\bm{\beta}}_{p}$ satisfies $\bm{E}^{\prime}(\theta)\hat{\bm{\beta}}_{p}=\bm{P}(\theta)\hat{\tilde{\bm{h}}}_{p}^{(m,m)}$, where

[TABLE]

Then, we can convert the function to log-likelihood to obtain

[TABLE]

where $\bm{M}(\theta)\triangleq(\bm{I}_{D_{m}}-\bm{P}(\theta))^{H}(\bm{I}_{D_{m}}-\bm{P}(\theta))$, and further $\bm{M}(\theta)=\bm{I}_{D_{m}}-\bm{P}(\theta)$. The estimator in (64) can also be implemented as in (35) utilizing the properties of the trace operation, which reduces the computational complexity since the second matrix is calculated only once during the search.

Appendix B Covariance Matrix in SEKF

The observation $\bm{f}_{k}^{(m)}$ in (36) can also be written as

[TABLE]

where we drop the superscripts for the sake of simplicity. Its mean is $\mathbb{E}\{\bm{f}_{k}\}=\text{vec}\{\bm{R}_{f}(\theta_{k})\}$, and covariance matrix $\bm{Q}_{k}$ is

[TABLE]

While the second term is known, the first term should be calculated. Using (65), it can be written as

[TABLE]

where $\bm{X}(p_{1},p_{2})$ can be written in two equivalent forms:

[TABLE]

Note that each entry of $\bm{X}(p_{1},p_{2})$ is a product of four variables, without summations. Since the channel is complex Gaussian, the identity $\mathbb{E}\left\{aa^{*}bb^{*}\right\}=\mathbb{E}\left\{aa^{*}\right\}\mathbb{E}\left\{bb^{*}\right\}+\mathbb{E}\left\{ab^{*}\right\}\mathbb{E}\left\{a^{*}b\right\}$ can be used, where $a$ and $b$ are complex Gaussian random variables. The essence of this identity is that different groups are assumed as independent for each term, e.g., $a$ and $b^{*}$ are assumed as independent of $a^{*}$ and $b$ in the second term. The two equivalent forms of $\bm{X}(p_{1},p_{2})$ in (68) and (69) can be used for these two terms. Consequently,

[TABLE]

Using (70), (67), and (66), the expression in (38) is obtained.

Appendix C Analytical NMSE Expression for BA-LS

For the denominator in the NMSE expression in Section VII-B1, $\mathbb{E}\{||\tilde{\bm{h}}_{p}^{(m,m)}||_{2}^{2}\}=\text{tr}(\tilde{\bm{R}}_{p}^{(m,m)})$ by definition. For the numerator, we define the error vector as $\bm{e}\triangleq\hat{\tilde{\bm{h}}}_{p}^{(m,m)}-\tilde{\bm{h}}_{p}^{(m,m)}$, which is evident from (20). Its terms are independent and the total energy of $\bm{e}$ is equal to the sum of energies of its terms. We assume that all the symbols have unit power, and $(\bm{s}^{(m)})^{H}\bm{s}^{(m)}$ is deterministic and equal to $N_{\text{F}}$. In addition, $\bm{s}^{(m^{\prime})}$ are uncorrelated random data for $m^{\prime}\neq m$, which yields $\mathbb{E}\{|(\bm{s}^{(m)})^{H}\bm{s}^{(m^{\prime})}|^{2}\}=N_{\text{F}}$. Then, the energy $\mathbb{E}\{||\bm{e}||_{2}^{2}\}=\text{tr}(\mathbb{E}\{\bm{e}\bm{e}^{H}\})$ of the error is written as (71) using (20), and simplified as (72) by inspecting (12).

[TABLE]