Fast Accurate Beam and Channel Tracking for Two-dimensional Phased Antenna Arrays

TL;DR

This paper develops fast, accurate 2D beam and channel tracking algorithms for mmWave systems with large phased arrays, reducing overhead and approaching theoretical limits, especially for dynamic channels.

Contribution

It introduces novel tracking algorithms optimized for different channel dynamics that minimize exploration overhead and achieve near-optimal accuracy.

Findings

Algorithms approach the Cramer-Rao lower bound for certain channels.

Simulation shows improved tracking accuracy and speed.

Efficient algorithms with low computational complexity.

Abstract

The sparsity and the severe attenuation of millimeter-wave (mmWave) channel imply that highly directional communication is needed. The narrow beam produced by large array requires accurate alignment, which is difficult to achieve when serving fast-moving users. In this paper, we focus on accurate two-dimensional (2D) beam and channel tracking problem aiming at minimizing exploration overhead and tracking error. Using a typical frame structure with periodic exploration and communication, a proven minimum overhead of exploration is provided first. Then tracking algorithms are designed for three types of channels with different dynamic properties. It is proved that the algorithms for quasi-static channels and channels in Dynamic Case I are optimal in approaching the minimum Cramer-Rao lower bound (CRLB). The computational complexity of our algorithms is analyzed showing their efficiency, and simulation results verify their advantages in both tracking error and tracking speed.

Tables

Table 1. TABLE I: Summary of the main acronyms.

Cramér-Rao lower bound	CRLB
analog beamforming	ABF
exploring beamforming vector	EBV
exploration and communication cycle	ECC
direction parameter vector	DPV
exploring beamforming matrix	EBM
joint beam and channel tracking	JBCT
recursive beam tracking	RBT

Table 2. TABLE II: The asymptotically optimal set of exploration offsets in Quasi-static Case.

${\widetilde{𝚫}}_{S,1}^{\ast }$	${\widetilde{𝚫}}_{S,2}^{\ast }$	${\widetilde{𝚫}}_{S,3}^{\ast }$
${[-0.0963,0.5098]}^{\text{T}}$	${[-0.2906,-0.2906]}^{\text{T}}$	${[0.5098,-0.0963]}^{\text{T}}$

Table 3. TABLE III: The asymptotically optimal set of exploration offsets in Dynamic Case I.

${\widetilde{𝚫}}_{DI,1}^{\ast }$	${\widetilde{𝚫}}_{DI,2}^{\ast }$	${\widetilde{𝚫}}_{DI,3}^{\ast }$
${[0.5486,0.2451]}^{\text{T}}$	${[-0.5462,0.2482]}^{\text{T}}$	${[-0.0012,-0.6837]}^{\text{T}}$

Table 4. TABLE IV: Number of required complex operations in each ECC.

3GPP NR	Algorithm 2	Algorithm 1 and Algorithm 3	Extended Kalman filter	Compressed sensing
6	28	45	1427	129088

Equations

$${\text{h}}_k=\eta \left({\theta }_k,{\varphi }_k\right){\beta }_k^c\text{a}({\text{x}}_k),$$

$$\text{a}({\text{x}}_k)={\text{a}}_1\left(x_{k,1}\right)\otimes {\text{a}}_2\left(x_{k,2}\right)$$

$${\text{a}}_1(x_{k,1})\triangleq {\left[1,e^{j2\pi \frac{x_{k,1}}{M}},\cdots ,e^{j2\pi \frac{M-1}{M}{x}_{k,1}}\right]}^{\text{T}}$$

$${\text{a}}_2(x_{k,2})\triangleq {\left[1,e^{j2\pi \frac{x_{k,2}}{N}},\cdots ,e^{j2\pi \frac{N-1}{N}{x}_{k,2}}\right]}^{\text{T}},$$

$$\beta \left({\text{x}}_k\right)\triangleq \eta \left({\text{x}}_k\right){\beta }_k^c,$$

$${\text{h}}_k=\beta \left({\text{x}}_k\right)\text{a}({\text{x}}_k).$$

$${\mathbf{\nu }}_{k,i}=\beta \left({\text{x}}_k\right){\text{w}}_{k,i}^{\text{H}}\text{a}({\text{x}}_k)\text{s}+{\mathbf{\zeta }}_{k,i}.$$

$$y_{k,i}={\mathbf{\nu }}_{k,i}\frac{{\text{s}}^{\text{H}}}{|\text{s}|}=$$

$$=$$

$${\text{y}}_k=|\text{s}|\beta \left({\text{x}}_k\right){\text{W}}_k^{\text{H}}\text{a}({\text{x}}_k)+{\text{z}}_k.$$

$$\mathcal{B}\left({\text{x}}_t\right)\triangleq \left(x_{t,1}-1,x_{t,1}+1\right)\times \left(x_{t,2}-1,x_{t,2}+1\right).$$

$${\text{W}}_k=$$

$${\hat{\mathbf{\psi }}}_k=$$

$$\underset{{\Xi }_k}{\min }$$

$$\mathbb{E}\left[{\hat{\text{h}}}_k\right]={\text{h}}_k,$$

$$\text{eqref}e{q}_{o}bservatio{n}_{v}ector,\text{eqref}e{q}_{F}c,\text{eqref}e{q}_{F}e,\text{vspace}-5mm$$

$${\text{w}}_{k,i}=\frac{1}{\sqrt{MN}}\text{a}\left({\mathbf{\omega }}_{k,i}\right),$$

$${\mathbf{\omega }}_{k,i}={\hat{\text{x}}}_{k-1}+{\mathbf{\Delta }}_i,i=1,\cdots ,q.\text{vspace}-0mm$$

$${\text{w}}_{k,i}=\frac{1}{\sqrt{MN}}\text{a}\left({\hat{\text{x}}}_{k-1}+{\mathbf{\Delta }}_i\right),i=1,\cdots ,q.$$

$$\underset{\Xi }{\min }$$

$$\text{eqref}e{q}_{o}bservatio{n}_{v}ector,\text{eqref}e{q}_{F}c,\text{eqref}e{q}_{F}e,\text{eqref}e{q}_{c}onstrant1,\text{eqref}e{q}_{s}v.$$

$$p_S({\text{y}}_k|\mathbf{\psi },{\text{W}}_k)=\frac{1}{{\pi }^3{\sigma }_z^6}{e}^{-\frac{{\Vert {\text{y}}_k-|\text{s}|\beta {\text{W}}_k^{\text{H}}\text{a}(\text{x})\Vert }_2^2}{{\sigma }_z^2}}.$$

$$\frac{1}{MN}\mathbb{E}\left[{\Vert {\hat{\text{h}}}_k-\text{h}\Vert }_2^2\right]$$

$$\ge$$

$$\triangleq$$

$$\text{V}\triangleq \frac{\partial \text{h}}{\partial {\mathbf{\psi }}^{\text{T}}}=$$

$$=$$

$${\text{I}}_S(\mathbf{\psi },{\text{W}}_l)$$

$$=\frac{2{|\text{s}|}^2}{{\sigma }_z^2}\text{Re}\left\{{\text{V}}^{\text{H}}{\text{W}}_l{\text{W}}_l^{\text{H}}\text{V}\right\}.$$

$$\underset{k\to +\infty }{\lim }{\hat{\text{x}}}_k=\text{x}.\text{vspace}-0mm$$

$$\underset{k\to +\infty }{\lim }{\text{W}}_k=\text{W}={\left[{\text{w}}_1,{\text{w}}_2,{\text{w}}_3\right]}^{\text{T}},\text{vspace}-0mm$$

$${\text{w}}_i\triangleq \frac{1}{\sqrt{MN}}\text{a}\left(\text{x}+{\mathbf{\Delta }}_{S,i}\right),i=1,2,3\text{vspace}-0mm$$

$$\underset{k\to +\infty }{\lim }k{C}_S^t(\mathbf{\psi },{\text{W}}_1,\cdots ,{\text{W}}_k)$$

$$=$$

$$C_S^{\min }(\mathbf{\psi })=$$

Full text

I Introduction

Millimeter-wave (mmWave) mobile communication is currently a hot topic due to its much wider bandwidth compared with the sub-6GHz spectrum. In mmWave channels, the much higher frequency leads to severe propagation loss, atmospheric absorption, penetration loss and other obstructions [1]. Fortunately, the shorter wavelength in the mmWave band allows the deployment of a larger antenna array, providing a considerable beamforming gain to compensate for the path loss [2, 3, 4, 5, 6]. For a hybrid or analog beamforming (ABF) system as cost-efficient ways to obtain this array gain, misalignment of beam direction may not only degrade the effective receiving power, increasing mutual interference, but also lead to the loss of beam observations due to the users’ mobility, especially in fast-varying environments [7, 8]. Therefore, accurate beam tracking is crucial for serving fast-moving users in mmWave mobile communication system.

In this paper, we will focus on the problem in ABF. Since only one RF chain connected with the antennas via programmable phase shifters is available in ABF, only one set of phase shifts can be applied (forming a so-called exploring beamforming vector (EBV) in this paper) and one dimension of the multiple-antenna channel can be observed at a certain time. Hence, in order to estimate the direction and the gain of the beam, the transceiver needs to try several different EBVs one by one. These EBVs can have a significant impact on tracking performance [9, 10, 11].

Although there already exist some beam tracking methods in [12, 13, 14, 15, 7, 16], which utilized historical exploring directions and observations to obtain current estimates, the EBVs were not optimized in those tracking algorithms. While beamforming resulting in the highest combining signal-to-noise ratio (SNR) is the best for data transmission, it is not the best for tracking accuracy [10, 11]. Optimal design of EBVs is necessary to achieve as accurate beam alignment as possible in mmWave mobile communication.

In [9], a beam tracking algorithm was proposed, trying to optimize the EBVs, assuming that the channel gain is known. In [10], the authors started to jointly track the channel gain and the beam direction with optimal EBVs. In [11], the optimization of the EBVs is converted to a convex problem and solved by online optimization toolboxes. Despite the progress, only one-dimensional (1D) array is supported for beam tracking optimization in these works. However, in most mobile applications, two-dimensional (2D) arrays are necessary, not only for providing much higher array gain, but also for supporting both horizontal and vertical beam direction variation [17, 18]. It brings huge challenges when extending the optimal EBVs in [10, 11] to 2D arrays. The optimization of the exploring directions in [10] relies on the symmetrical property of the two exploring beams, however, keeping this symmetrical property may result in more overhead when extended to 2D tracking. As for the algorithm in [11], the objective function versus the EBVs becomes extremely complicated and quite different in 2D tracking compared with 1D tracking, leading to the failure of the previously used optimization method in 1D system. To the best of the authors’ knowledge, there is no work on the design of optimal EBVs for 2D array yet.

In this paper, we focus on the design of the optimal EBVs and the accurate single-path tracking algorithms based on 2D phased antenna array. The widely used frame structure [10, 8, 19] is adopted here. As shown in Fig. 1, the transceiver periodically works in exploration and communication mode. In the exploration stage of each exploration and communication cycle (ECC), the transmitter sends a pre-defined pilot sequence for $q$ times. At each time, the receiver forms one exploring beam pointing in one direction to observe the channel. Then the channel gain and the direction of the incoming beam are estimated with the $q$ observations of the channel. In the communication stage of each ECC, the beam is aligned in the current estimated direction, and the current estimated channel gain will be used for the subsequent process. Based on this structure, the following questions are to be answered:

1. What is the minimum exploration overhead $q$ in each ECC for 2D tracking?

2. How to determine the $q$ exploring directions based on the channel observations in previous ECCs for 2D tracking?

3. How to track the 2D beam direction and the channel gain for different time-varying channels, e.g., from quasi-static channels to fast-fading channels?

4. How is the accuracy, convergence and stability of the tracking algorithm?

Following these questions, we summarize the main contributions of this paper as below:

1. Based on a reasonable EBV constraint, it is proved that the minimum exploration overhead counted by the number of exploring directions is $q=3$, for a unique solution of the 2D beam direction and the channel gain within only one ECC, while simple extension from 1D to 2D tracking will need $q=4$.

2. Dynamic beam and channel tracking strategies for three different time-varying channels (called Quasi-static Case, Dynamic Case I and Dynamic Case II in this paper) are proposed and optimized. The salient advantages of these tracking algorithms are given below:

i) In Quasi-static Case (channels with quasi-static beam direction and channel gain), the optimal exploration offsets are derived. Also, a joint beam direction and channel gain tracking algorithm is proposed, and the tracking error is proved to converge to the minimum Cramér-Rao lower bound (CRLB).

ii) In Dynamic Case I (channels with quasi-static beam direction and fast-fading channel gain), the Rayleigh fading channel is studied as a special case in this paper. The optimal exploration offsets are obtained and an algorithm for beam (only) tracking is proposed, which is proved to converge and achieve the minimum CRLB on the beam direction.

iii) In Dynamic Case II (channels with fast-changing beam direction and channel gain), a joint tracking algorithm of the beam direction and the channel gain is proposed with faster and more accurate performance.

3. The impact of the antenna pattern on the tracking algorithms and the performance is taken into account, showing that the proposed algorithms are suitable for practical implementations.

Part of this work was presented in our conference paper [20], while the main difference and novelty of this paper lies in the following four aspects: 1) tracking with a general antenna pattern rather than a simple isotropic pattern: we consider a more general direction-dependent antenna element pattern here rather than an isotropic pattern; 2) tracking for different types of time-varying channels: in addition to slow-fading channels, fast-fading channels are also studied in this paper and the corresponding tracking strategy is proposed and optimized, while the prior algorithm does not support fast-fading channels; 3) the complexity analysis: we analyze the computational complexity of the proposed tracking algorithms while this was missing in our previous work; 4) more rigorous conclusions and more complete proofs: we provide more rigorous lemmas and theorems with more complete proofs in this paper, correcting the corresponding flaws in the previous work.

The remaining part of this paper is organized as follows: the system model is described in Section II. In Section III, the tracking problem with some constraints is formulated. Then the minimum exploration overhead of joint 2D beam and channel tracking is given in theory in Section IV. In Section V and Section VI, the tracking problems for Quasi-static Case (Section V) and Dynamic Case I (Section VI) are studied separately. The tracking performance bounds are derived and corresponding tracking algorithms are developed with convergence and optimality analysis. In Section VII, a tracking algorithm is developed for Dynamic Case II. Then the complexity analysis of these algorithms is given in Section VIII. Section IX presents numerical results to verify the performance of our proposed algorithms.

Notations: We use lower case letters such as $a$ and a to denote scalars and column vectors. Respectively, $\lvert\textbf{a}\rvert$ and $\left\|\textbf{a}\right\|_{2}$ represent the modulus and 2-norm of the vector a. Upper case boldface letters, e.g., A, are used to denote matrices. The superscript $\bar{\left(\cdot\right)}$, $\left(\cdot\right)^{\text{T}}$, $\left(\cdot\right)^{\text{H}}$ are utilized to denote conjugate, transpose and conjugate-transpose. For a matrix A, its inverse, pseudo-inverse and determinant are written as $\textbf{A}^{-1}$, $\textbf{A}^{+}$ and $\lvert\textbf{A}\rvert$. The identity matrix of order $q$ is denoted by $\textbf{J}_{q}$. Let $\mathcal{CN}(\mu,\sigma^{2})$ represent the symmetric complex Gaussian distribution with mean $\mu$ and variance $\sigma^{2}$, and $\mathcal{N}(\mu,\sigma^{2})$ stand for the real Gaussian distribution with mean $\mu$ and variance $\sigma^{2}$. The Kronecker product is represented as $\otimes$. The statistical expectation is denoted by $\mathbb{E}\left[\cdot\right]$. The real (imaginary) part is represented as $\text{Re}\left\{\cdot\right\}$ $\left(\text{Im}\left\{\cdot\right\}\right)$. The natural logarithm of a scalar $y$ is obtained by $\log\left(\cdot\right)$ and the phase angle of a complex number $z$ is written as $\angle z$. The main acronyms used in this paper are summarized in TABLE I.

II System Model

II-A System Configuration

We consider a mmWave receiver111 Note that tracking is needed at both the transmitter and the receiver. However, considering the transmitter-receiver reciprocity, the tracking of both sides have similar designs. Hence, we focus on beam and channel tracking at the receiver side. equipped with a planar phased antenna array, as shown in Fig. 2. The planar array consists of $M\times N$ antenna elements that are placed in a rectangular area, where $M$($N$) antenna elements are evenly distributed along $x$-axis ($z$-axis) with a distance $d_{1}$ ($d_{2}$) between neighboring elements. These antenna elements are connected to the same RF chain via programmable phase shifters.

Single RF chain of ABF makes a constraint that only one beam can be formed at any time and hence the receiver has to work alternatively in exploration and communication mode, resulting in a frame structure of periodic ECC. The angle of arrival (AoA) and the channel gain are assumed to be constant in each ECC and may change in different ECCs. In the exploration stage of one ECC, the transmitter sends a pre-defined pilot sequence s for $q$ times, where $\textbf{s}=\left[s_{1},\cdots,s_{L_{s}}\right]\in\mathbb{C}^{1\times L_{s}}$ contains $L_{s}$ same symbols. At each time, the receiver forms one exploring beam pointing in one direction to observe the channel. Then the channel gain and the direction of the incoming beam are estimated according to the $q$ observations obtained in the current and previous ECCs. In the communication stage of each ECC, the beam is aligned in current estimated direction, and the current estimated channel gain will be used for the subsequent process.

II-B Channel Model

In mmWave outdoor communication, the scattering is not rich and the number of effective propagation paths is usually limited [1, 21]. Besides, the beam formed by a large array in the mmWave system is quite narrow and the interaction between multi-path is relatively weak [22]. In other words, the incoming paths are usually sparse in space, making it possible to track each path independently. Hence, we focus on the method for tracking one path, while different paths can be tracked separately by using the same method.

In $k$-th ECC, the direction of the incoming beam is denoted by ($\theta_{k},\phi_{k}$), where $\theta_{k}\in[-\frac{\pi}{2},\frac{\pi}{2})$ is the elevation AoA and $\phi_{k}\in[0,\pi)$ is the azimuth AoA. Then the channel vector of this path during $k$-th ECC is

[TABLE]

where $\eta\left(\theta,\phi\right)$ is the direction-dependent antenna gain (antenna pattern) of each element,222 Note that the antenna patterns of different elements in the array may not be exactly the same and need to be carefully calibrated. However, considering that it is not the focus of this paper, we assume the same patterns for all the elements here. $\beta_{k}^{c}$ is the complex channel gain, \textbf{x}_{k}\triangleq\left[x_{k,1},x_{k,2}\right]^{\text{T}}=\big{[}\frac{Md_{1}\cos(\theta_{k})\cos(\phi_{k})}{\lambda}, \frac{Nd_{2}\sin(\theta_{k})}{\lambda}\big{]}^{\text{T}}is the direction parameter vector (DPV) determined by ($\theta_{k},\phi_{k}$),

[TABLE]

is the 2D steering vector with

[TABLE]

and $\lambda$ is the wavelength. For the convenience of expression, the antenna gain $\eta\left(\theta_{k},\phi_{k}\right)$ is denoted by $\eta\left(\textbf{x}_{k}\right)$ hereinafter.

Define the equivalent channel gain in $k$-th ECC as below:

[TABLE]

then the channel vector in (1) can be rewritten as

[TABLE]

II-C RF and Baseband Preprocessing

Synchronization in both carrier frequency and symbol timing is necessary in mmWave wireless communications. In the initial beam estimation stage in Fig. 1, the carrier frequency synchronization information can be obtained and estimated, the residual error of which can be converted to the time-varying phase of the equivalent channel gain in (5). As for the symbol timing, since it changes much slower, it can be estimated and tracked much more easily both in the initial beam estimation stage and the tracking stage. There already exists a lot of work on the synchronization algorithms in mmWave communication systems [23, 24, 25]. To make the research goals more focused, we assume perfect synchronization in this paper. Future work may be needed to further study the impact of residual synchronization error on the beam tracking performance.

Next, we will focus on the receiving beamforming based on the perfect synchronization assumption above. Let $\textbf{w}_{k,i}\in\mathbb{C}^{MN\times 1}$ be the EBV for receiving the pilot sequence s the $i$-th $(i=1,\cdots,q)$ time in $k$-th ECC. The entries of $\textbf{w}_{k,i}$ are of the same amplitude with \bigg{|}\left[\textbf{w}_{k,i}\right]_{l}\bigg{|}=\frac{1}{\sqrt{MN}}, where $\left[\textbf{w}_{k,i}\right]_{l}$ denotes the $l$-th element of $\textbf{w}_{k,i}$. After phase shifting and combining, the $i$-th received sequence in $k$-th ECC at the baseband output of the RF chain is given by

[TABLE]

where $\boldsymbol{\zeta}_{k,i}\in\mathbb{C}^{1\times L_{s}}$ is the receiving noise vector.

By match filtering on the sequence $\boldsymbol{\nu}_{k,i}$, the $i$-th observation in $k$-th ECC is given below:

[TABLE]

where $z_{k,i}\triangleq\boldsymbol{\zeta}_{k,i}\frac{\textbf{s}^{\text{H}}}{\lvert\textbf{s}\rvert}$ is an additive noise, which is modeled as i.i.d. Gaussian distributed in this paper, i.e., $z_{k,i}\sim\mathcal{CN}\left(0,\sigma_{z}^{2}\right)$. This assumption is certainly held when the receiving noise vector $\boldsymbol{\zeta}_{k,i}$ is i.i.d. Gaussian distributed. Besides, even when the noise vector $\boldsymbol{\zeta}_{k,i}$ is non-Gaussian, if the real and imaginary parts of the elements in $\boldsymbol{\zeta}_{k,i}$ are i.i.d, the observation noise $z_{k,i}$ can also be regarded as i.i.d. Gaussian distributed as long as the pilot sequence length $L_{s}$ is sufficiently large, according to the central limit theorem [26].

Let $\textbf{W}_{k}\triangleq\left[\textbf{w}_{k,1},\ldots,\textbf{w}_{k,q}\right]$, $\textbf{z}_{k}\triangleq\left[z_{k,1},\ldots,z_{k,q}\right]^{\text{T}}$ and $\textbf{y}_{k}\triangleq\left[y_{k,1},\ldots,y_{k,q}\right]^{\text{T}}$ denote the exploring beamforming matrix (EBM), the noise vector and the observation vector respectively. Then we can rewrite (8) as follows:

[TABLE]

II-D Tracking Loop

As shown in Fig. 1, an initial estimate $\hat{\beta}_{0}=\hat{\beta}_{0}^{\text{re}}+j\hat{\beta}_{0}^{\text{im}}$ and $\hat{\textbf{x}}_{0}=\left[\hat{x}_{0,1},\hat{x}_{0,2}\right]^{\text{T}}$ can be obtained in the beam estimation stage. It is assumed in this paper that the initial beam estimator can output an estimate $\hat{\textbf{x}}_{0}$ falling within the main lobe of $\textbf{x}_{0}$, i.e., $\hat{\textbf{x}}_{0}\in\mathcal{B}\left(\textbf{x}_{0}\right)$, where $\mathcal{B}\left(\textbf{x}_{t}\right)$ denotes the main lobe of an arbitrary DPV \textbf{x}_{t}=\big{[}x_{t,1},x_{t,2}\big{]}^{\text{T}}, given by

[TABLE]

Then our tracking starts from this initial estimate $\hat{\textbf{x}}_{0}$ to find more accurate beam directions. It is worth pointing out that the main lobe in x domain in (10) has been normalized to a square with twice the unit length of each side and centered at the DPV $\textbf{x}_{t}$ after the transformation from the angle domain to x domain. Hence, the main lobe size in x domain remains unchanged even if the antenna size $M,\,N$ scale.

In the exploration stage of $k$-th ECC, the receiver needs to choose an EBM $\textbf{W}_{k}$ based on historical observation vectors $\textbf{y}_{1},\cdots,\textbf{y}_{k-1}$ along with the corresponding EBMs $\textbf{W}_{1},\cdots,\textbf{W}_{k-1}$. The new observation $\textbf{y}_{k}$ can be obtained by applying $\textbf{W}_{k}$. Then the estimate $\hat{\boldsymbol{\psi}}_{k}\triangleq\left[\hat{\beta}_{k}^{\text{re}},\hat{\beta}_{k}^{\text{im}},\hat{x}_{k,1},\hat{x}_{k,2}\right]^{\text{T}}$ of the channel parameter vector \boldsymbol{\psi}_{k}\triangleq\big{[}\text{Re}\left\{\beta\left(\textbf{x}_{k}\right)\right\},\text{Im}\left\{\beta\left(\textbf{x}_{k}\right)\right\}, x_{k,1},{x}_{k,2}\big{]}^{\text{T}} is obtained by using all observation vectors available and the corresponding EBMs. The whole tracking loop is given in Procedure 1 and the focus of this paper lies in Step 3 and Step 6.

From a control system perspective, ${\boldsymbol{\psi}}_{k}$ is the system state, $\hat{{\boldsymbol{\psi}}}_{k}$ is the estimate of the system state, the EBM $\textbf{W}_{k}$ is the control action and $\textbf{y}_{k}$ is a noisy observation non-linearly determined by the system state and the control action. Hence, the task of a tracking design is to find the following strategy:

[TABLE]

where $\textbf{F}_{k}^{c}$ denotes the control function and $\textbf{F}_{k}^{e}$ denotes the estimation function in $k$-th ECC.

III Problem Formulation

Let $\Xi_{k}=\left\{\textbf{F}_{k}^{c},\textbf{F}_{k}^{e}\right\}$ denote the set of beam and channel tracking schemes in $k$-th ECC. Then the optimal beam and channel tracking problem minimizing the mean square error (MSE) of the channel vector estimate is formulated as:

[TABLE]

where the constraint (14) ensures that $\hat{\textbf{h}}_{k}\triangleq\hat{\beta}_{k}\textbf{a}\left(\hat{\textbf{x}}_{k}\right)$ is an unbiased estimate of the channel vector $\textbf{h}_{k}=\beta\left(\textbf{x}_{k}\right)\textbf{a}\left(\textbf{x}_{k}\right)$. It is worth explaining the following two points. First, an unbiased estimator may not be the best estimator that achieves the minimum MSE. Nevertheless, such an optimal estimator with no constraints is hard to obtain and hence we add this unbiasedness constraint. Second, we only need to guarantee the unbiasedness of $\hat{\textbf{h}}_{k}$, as the objective function in (13) is the MSE of the channel vector. The estimate of the equivalent channel gain and the DPV, i.e., $\hat{\beta}_{k}$ and $\hat{\textbf{x}}_{k}$, can be biased.

Problem (13) is challenging to be solved optimally due to the following reasons:

1. It is a partially observed Markov decision process (POMDP) which generally has not been solved optimally [27, 28].

2. There are $M\times N$ phase shifts to adjust in each EBV $\textbf{w}_{k,i}$. This makes the optimization of the EBV too complicated due to the joint design of so many phase shifts, especially when the array size $M\times N$ grows large.

3. To obtain $\hat{\boldsymbol{\psi}}_{k}$ in $k$-th ECC, $k$ EBMs, i.e., $\textbf{W}_{1},\cdots,\textbf{W}_{k}$, need to be designed, making it difficult to optimize so many beamforming matrices simultaneously as $k$ increases.

4. The time-varying features of the channel vector in (6) restrict the tracking algorithm and the system performance. Thus, it is hard to design an optimal tracking method for a general channel model.

These challenges above make it extremely difficult to solve this problem optimally. Hence, we add some reasonable constraints in this paper to take the first step of the optimal tracking policy:

III-A The EBV constraint

As it is complicated to obtain the optimal $M\times N$ phase shifts in general for each EBV, we use steering vectors to design the EBVs,

[TABLE]

where $\boldsymbol{\omega}_{k,i}\triangleq\left[{\omega}_{k,i1},{\omega}_{k,i2}\right]^{\text{T}}$ denotes the $i$-th exploring direction vector in $k$-th ECC. This ensures that only two variables need to be designed for each EBV.

III-B The exploring direction constraint

Although the exploring direction vector $\boldsymbol{\omega}_{k,i}$ in (15) can be of any form, however, considering the tracking accuracy, it is better to make sure that $\boldsymbol{\omega}_{k,i}$ falls within the main lobe of the DPV $\textbf{x}_{k}$ in (10). Thus, it is reasonable to choose exploring directions near the recently estimated direction $\hat{\textbf{x}}_{k-1}$. For this purpose, we use such an architecture in this paper. That is, the $i$-th exploring direction vector in $k$-th ECC, i.e., $\boldsymbol{\omega}_{k,i}$, is determined by the previous estimate of the DPV plus an exploration offset $\boldsymbol{\Delta}_{k,i}$. Considering the design of the offsets that change in different ECCs is also very complicated, we adopt fixed exploration offsets $\boldsymbol{\Delta}_{i}(i=1,\cdots,q)$ in this paper:

[TABLE]

Therefore, the EBV in (15) can be rewritten as

[TABLE]

III-C The time-varying channel constraint

The time-varying channel vector in (1) is determined by three parts: the antenna gain $\eta\left(\textbf{x}_{k}\right)$, the channel gain $\beta_{k}^{c}$ and the DPV $\textbf{x}_{k}$. Since the change of the antenna gain $\eta\left(\textbf{x}_{k}\right)$ depends on the DPV $\textbf{x}_{k}$ for a given antenna element pattern, we only consider the change of the DPV $\textbf{x}_{k}$ and the channel gain $\beta_{k}^{c}$ when exploring the properties of the time-varying channels. As the user motion characteristics can be quite different in various situations [29, 30, 31], both of the DPV $\textbf{x}_{k}$ and the channel gain $\beta_{k}^{c}$ may change slowly or fast. Therefore, four possible cases exist, which correspond to four different practical scenarios and can be modeled as follows:

• •

Quasi-static Case: $\textbf{x}_{k}\approx\textbf{x},\beta_{k}^{c}\approx\beta^{c}$

When both $\textbf{x}_{k}$ and $\beta_{k}^{c}$ change slowly, e.g., the user keeps static or quasi-static in a room, the antenna gain $\eta\left(\textbf{x}_{k}\right)$ and the equivalent channel gain $\beta\left(\textbf{x}_{k}\right)$ defined in (5) also change slowly. The channel in this case can be seen as approximately fixed. For the sake of convenience, we assume that $\beta\left(\textbf{x}_{k}\right)=\beta=\beta^{\text{re}}+j\beta^{\text{im}}$, $\textbf{x}_{k}=\textbf{x}=\left[x_{1},x_{2}\right]^{\text{T}}$ in this case.

• •

Dynamic Case: $\textbf{x}_{k}\approx\textbf{x},\beta_{k+1}^{c}\neq\beta_{k}^{c}$

For channels that $\textbf{x}_{k}$ changes slowly while $\beta_{k}^{c}$ changes fast, e.g., a person walks at a fast pace in a room, the beam direction can be seen as approximately fixed, i.e., $\textbf{x}_{k}=\textbf{x}$[29, 30]. To distinguish from other dynamic scenarios, this case is called Dynamic Case I.

• •

Dynamic Case: $\textbf{x}_{k+1}\neq\textbf{x}_{k},\beta_{k}^{c}\approx\beta^{c}$

This case requires that the beam direction changes fast while the channel gain keeps static or quasi-static. However, in real mmWave channels, the fast change of the beam direction usually leads to the fast change of the channel gain since the propagation paths change. This case exists only when the user rotates around the base station (BS) exactly in a circle in line of sight (LOS) channels. This is not the usual case and not studied in this paper.

• •

Dynamic Case: $\textbf{x}_{k+1}\neq\textbf{x}_{k},\beta_{k+1}^{c}\neq\beta_{k}^{c}$

Both the beam direction and the channel gain in this case change fast, which happens in most fast-moving scenarios except Dynamic Case I, e.g., an unmanned aerial vehicle (UAV) flies in the sky [31]. To distinguish from Dynamic Case I, we call it Dynamic Case II.

It would be helpful to explain the following two aspects. First of all, extra algorithms need to be introduced to efficiently classify the channels according to the time-varying features. However, to make the research goals more focused in this paper, we leave the details of this classification in future work. Second, this paper only exploits the independent variation properties of $\beta_{k}^{c}$ and $\textbf{x}_{k}$ in these four cases to obtain theoretical results. While in real mmWave channels, the variation of the channel gain $\beta_{k}^{c}$ and the DPV $\textbf{x}_{k}$ might be interrelated [32], which are supposed to be jointly taken into account in future work.

With the above-mentioned EBV constraint, the exploring direction constraint and the time-varying channel constraint, the beam and channel tracking problem in (13) can be reformulated as:

[TABLE]

IV How Many Explorations Are Needed In each ECC?

Before delving into the detailed tracking process in (18), we will first study the number of explorations needed in this section.

To estimate ${\boldsymbol{\psi}}_{k}$, sufficient measurements from different exploring directions are required. For Quasi-static Case where ${\boldsymbol{\psi}}_{k}$ remains unchanged, i.e., ${\boldsymbol{\psi}}_{k}=\boldsymbol{\psi}\triangleq\left[\beta^{\text{re}},\beta^{\text{im}},x_{1},x_{2}\right]^{\text{T}}$, one exploration in each ECC is enough since sufficient measurements are available after quite a number of ECCs. Nevertheless, in dynamic case, only using one exploration in each ECC does not work well as ${\boldsymbol{\psi}}_{k}$ may change fast. Hence, it is necessary to ensure that the estimate can be obtained even by using the explorations in a single ECC. Then the question becomes: under the condition above, how many explorations are needed in each ECC?

With the constraint in (15), two explorations in each ECC are sufficient to jointly track the equivalent channel gain and the 1D beam direction according to [10]. When tracking the 2D direction, it is straight forward that four explorations are feasible by separately using two explorations to track each dimension of the 2D direction. However, using four explorations will lower the system efficiency since it will cost time resources for each exploration. Hence, we may ask that can we reduce the times of exploration, or what is the minimum number of explorations required?

Then the following lemma is proposed to help determine the minimum exploration overhead $q$ in each ECC:

Lemma 1**.**

If the EBVs are of the steering vector forms, i.e., $\textbf{w}_{k,i}=\frac{1}{\sqrt{MN}}\textbf{a}\!\left(\boldsymbol{\omega}_{k,i}\right)$, and the observation vector in (9) is noiseless, then

1) to obtain the unique solution of the channel parameter vector $\boldsymbol{\psi}_{k}$ within one ECC, the minimum exploration overhead is $q=3$ in each ECC;

2) to obtain the unique solution of the DPV $\textbf{x}_{k}$ within one ECC, the minimum exploration overhead is $q=3$ in each ECC.

Proof.

See Appendix A. ∎

Lemma 1 reveals that it is impossible to obtain the unique solution within one ECC when only using two explorations, whether we want to jointly estimate $\beta\left(\textbf{x}_{k}\right)$ and $\textbf{x}_{k}$ or just estimate $\textbf{x}_{k}$. If we use three explorations and design three appropriate exploring directions in each ECC, then the unique solution of the channel parameter vector $\boldsymbol{\psi}_{k}$ can be obtained. Hence, we set $q=3$ in this paper, i.e., the EBM $\textbf{W}_{k}=\left[\textbf{w}_{k,1},\textbf{w}_{k,2},\textbf{w}_{k,3}\right]$.

V Quasi-static Tracking: Performance Bound, Convergence and Optimality

In this section, we will focus on Quasi-static Case. As mentioned in Section III, in Quasi-static Case, $\boldsymbol{\psi}_{k}=\boldsymbol{\psi}=\left[\beta^{\text{re}},\beta^{\text{im}},x_{1},x_{2}\right]^{\text{T}}$ and $\textbf{h}_{k}=\textbf{h}\triangleq{\beta}\textbf{a}\left({\textbf{x}}\right)$. For a given channel parameter vector $\boldsymbol{\psi}$ and EBM $\textbf{W}_{k}$, the observation vector satisfies normal distribution with $\textbf{y}_{k}\sim\mathcal{CN}\left(\lvert\textbf{s}\rvert\beta\textbf{W}_{k}^{\text{H}}\textbf{a}(\textbf{x}),\sigma_{z}^{2}\textbf{J}_{3}\right)$. Hence, the conditional probability density function of $\textbf{y}_{k}$ is given by

[TABLE]

In this section, we will first provide the lower bound of the tracking error in Quasi-static Case. Then we develop a tracking algorithm and prove it can converge to the minimum CRLB with time.

V-A Cramér-Rao Lower Bound of The Tracking Error

The Cramér-Rao lower bound theory gives the lower bound of the unbiased estimation error [33]. Based on this, we introduce the following lemma to obtain the lower bound of the tracking error in Quasi-static Case:

Lemma 2**.**

In Quasi-static Case, given $\textbf{W}_{1},\cdots,\textbf{W}_{k}$, the MSE of the channel vector estimate in (18) is lower bounded as follows:

[TABLE]

where V is the Jacobian matrix given by

[TABLE]

and the Fisher information matrix $\textbf{I}_{S}(\boldsymbol{\psi},\textbf{W}_{l})$ is given by

[TABLE]

Proof.

See Appendix B. ∎

The CRLB in (20) is a function of the EBMs $\textbf{W}_{1},\ldots,\textbf{W}_{k}$. Since it is hard to optimize so many EBMs simultaneously, we will first try to find a lower bound of the CRLB under the constraint (17), and later design a tracking algorithm approaching this lower bound.

Consider any tracking algorithm under the constraint (17) that can converge to the DPV x, i.e.,

[TABLE]

Then the EBM $\textbf{W}_{k}$ also converges,

[TABLE]

where $\textbf{w}_{i}$ is given by

[TABLE]

with $\left\{\boldsymbol{\Delta}_{S,1},\boldsymbol{\Delta}_{S,2},\boldsymbol{\Delta}_{S,3}\right\}$ denoting the fixed set of exploration offsets in Quasi-static Case. Hence, the normalized CRLB (by multiplying $k$) converges as $k\to+\infty$:

[TABLE]

which is a function of $\boldsymbol{\psi},\,\textbf{W}$ and will be denoted as ${C}_{S}(\boldsymbol{\psi},\textbf{W})$.

According to (V-A), for a given channel (direction and gain), there exists an optimal EBM, which leads to the minimum normalized CRLB as a function of the channel parameter vector $\boldsymbol{\psi}$:

[TABLE]

Solving problem (27) yields $\textbf{W}_{S}^{*}=\left[\textbf{w}_{S,1}^{*},\textbf{w}_{S,2}^{*},\textbf{w}_{S,3}^{*}\right]$ with

[TABLE]

where $\left\{\boldsymbol{\Delta}_{S,1}^{*},\boldsymbol{\Delta}_{S,2}^{*},\boldsymbol{\Delta}_{S,3}^{*}\right\}$ denotes the optimal set of exploration offsets for a given array size and a given $\boldsymbol{\psi}$.

V-B Asymptotically Optimal Set of Exploration Offsets

In general, the minimum CRLB in (27) is a function of a set of system parameters including the equivalent channel gain $\beta$, the DPV x and the array size $M,\,N$. Hence, the optimal set of 2D exploration offsets should also be a function of these parameters. Since it is very hard to obtain the expression of this optimal set, we adopt numerical search to deal with this issue. However, as many parameters in (27) may affect the optimal result, numerical search has to be reconducted for different parameter sets, resulting in high complexity.

Fortunately, through our investigation, some useful properties of the minimum CRLB and the optimal set of exploration offsets are given to simplify the numerical search, as described in the following lemma:

Lemma 3**.**

In Quasi-static Case, the minimum CRLB ${C}_{S}^{\min}(\boldsymbol{\psi})$ and the optimal set of exploration offsets $\left\{\boldsymbol{\Delta}_{S,1}^{*},\boldsymbol{\Delta}_{S,2}^{*},\boldsymbol{\Delta}_{S,3}^{*}\right\}$ have the following three properties:

1) ${C}_{S}^{\min}(\boldsymbol{\psi})$, $\left\{\boldsymbol{\Delta}_{S,1}^{*},\boldsymbol{\Delta}_{S,2}^{*},\boldsymbol{\Delta}_{S,3}^{*}\right\}$ are invariant to the equivalent channel gain $\beta$;

2) ${C}_{S}^{\min}(\boldsymbol{\psi})$, $\left\{\boldsymbol{\Delta}_{S,1}^{*},\boldsymbol{\Delta}_{S,2}^{*},\boldsymbol{\Delta}_{S,3}^{*}\right\}$ are invariant to the DPV x;

3) ${C}_{S}^{\min}(\boldsymbol{\psi})$ converges as $\emph{M},\,\emph{N}\to+\infty$ and there exists a fixed set of exploration offsets that are unrelated to the array size $M,N$, denoted as $\left\{\widetilde{\boldsymbol{\Delta}}_{S,1}^{*},\,\widetilde{\boldsymbol{\Delta}}_{S,2}^{*},\widetilde{\boldsymbol{\Delta}}_{S,3}^{*}\right\}$, such that

[TABLE]

where $\widetilde{\textbf{W}}_{S}^{*}=[\tilde{\textbf{w}}_{S,1}^{*},\tilde{\textbf{w}}_{S,2}^{*},\tilde{\textbf{w}}_{S,3}^{*}]$ is obtained with

[TABLE]

Proof.

See Appendix C. ∎

Lemma 3 reveals that $\left\{\boldsymbol{\Delta}_{S,1}^{*},\boldsymbol{\Delta}_{S,2}^{*},\boldsymbol{\Delta}_{S,3}^{*}\right\}$ is only related to the array size $M,\,N$. Hence, the numerical search times can be reduced to one for a particular array size $M,\,N$. Numerically, we find later that even if $\left\{\boldsymbol{\Delta}_{S,1}^{*},\boldsymbol{\Delta}_{S,2}^{*},\boldsymbol{\Delta}_{S,3}^{*}\right\}$ may change for different array sizes, $\left\{\widetilde{\boldsymbol{\Delta}}_{S,1}^{*},\,\widetilde{\boldsymbol{\Delta}}_{S,2}^{*},\widetilde{\boldsymbol{\Delta}}_{S,3}^{*}\right\}$ can be used to take the place of $\left\{\boldsymbol{\Delta}_{S,1}^{*},\boldsymbol{\Delta}_{S,2}^{*},\boldsymbol{\Delta}_{S,3}^{*}\right\}$ as long as $M$ and $N$ are sufficiently large. Therefore, the numerical search times is reduced to one in the end. As $\left\{\widetilde{\boldsymbol{\Delta}}_{S,1}^{*},\,\widetilde{\boldsymbol{\Delta}}_{S,2}^{*},\widetilde{\boldsymbol{\Delta}}_{S,3}^{*}\right\}$ can be used to achieve the minimum CRLB when $M,N\to+\infty$, it is called the asymptotically optimal set of exploration offsets in Quasi-static Case in this paper.

By numerical search in the main lobe in (10), we can obtain one asymptotically optimal set of exploration offsets $\left\{\widetilde{\boldsymbol{\Delta}}_{S,1}^{*},\,\widetilde{\boldsymbol{\Delta}}_{S,2}^{*},\widetilde{\boldsymbol{\Delta}}_{S,3}^{*}\right\}$ in TABLE II and Fig. 3. It can be seen that the three exploring direction vectors do not form a regular triangle as the radiation pattern produced by (25) is not isotropic from different angles.

With this set in TABLE II, a general way to generate the EBM $\widetilde{\textbf{W}}_{S}^{*}$ is obtained by (29).

The set of exploration offsets $\left\{\widetilde{\boldsymbol{\Delta}}_{S,1}^{*},\,\widetilde{\boldsymbol{\Delta}}_{S,2}^{*},\widetilde{\boldsymbol{\Delta}}_{S,3}^{*}\right\}$ may become sub-optimal when the antenna size $M\times N$ is finite. To evaluate the robustness of this set of exploration offsets to finite array size, we adopt $\left\{\widetilde{\boldsymbol{\Delta}}_{S,1}^{*},\,\widetilde{\boldsymbol{\Delta}}_{S,2}^{*},\widetilde{\boldsymbol{\Delta}}_{S,3}^{*}\right\}$ to antenna arrays of limited size and compare the minimum CRLB with the CRLB achieved by $\left\{\widetilde{\boldsymbol{\Delta}}_{S,1}^{*},\,\widetilde{\boldsymbol{\Delta}}_{S,2}^{*},\widetilde{\boldsymbol{\Delta}}_{S,3}^{*}\right\}$ in TABLE II. As illustrated in Fig. 4, when the antenna number $M=N\geq 8$, we can approach the minimum CRLB with a relative error less than $0.1\%$ by using $\left\{\widetilde{\boldsymbol{\Delta}}_{S,1}^{*},\,\widetilde{\boldsymbol{\Delta}}_{S,2}^{*},\widetilde{\boldsymbol{\Delta}}_{S,3}^{*}\right\}$.

As a conclusion, it is practical to apply this asymptotically optimal set of exploration offsets to any antenna array with $M=N\geq 8$, any channel gain and any direction.

V-C Joint Beam and Channel Tracking

In the above subsections, we have provided a low-complexity numerical method to design the optimal exploration offsets and obtain the minimum CRLB, given that the DPV x is known. However, in a real tracking problem, the DPV x is unknown and the EBMs need to be adjusted dynamically. In addition, a sequence of optimal beamforming matrices only tells us what the minimum CRLB is, but it cannot tell us which tracking algorithm can achieve the minimum CRLB. In this subsection, we propose a specific tracking algorithm to approach the minimum CRLB.

The proposed tracker is motivated by the following maximum likelihood problem:

[TABLE]

Since $\textbf{y}_{1},\cdots,\textbf{y}_{k}$ are independently observed vectors, we can convert (30) as follows:

[TABLE]

This problem is somewhat similar to a kind of estimation problem with control in [34, Section 10.2], where a two-layer recursive algorithm is proposed and proved to converge to the minimum CRLB under several requirements. The difference is that the observation vector in [34, Section 10.2] is given by $\textbf{y}=\textbf{f}^{o}\left(\boldsymbol{\psi},\textbf{W}\right)+\textbf{z}$, where z is an i.i.d. Gaussian noise vector and $\textbf{f}^{o}\left(\boldsymbol{\psi},\textbf{W}\right)$ is a convex function with respect to $\boldsymbol{\psi}$. However, in our problem, $\textbf{f}^{o}\left(\boldsymbol{\psi},\textbf{W}\right)=\lvert s\rvert\beta\textbf{W}^{\text{H}}\textbf{a}\left(\textbf{x}\right)$ is non-convex, leading to the failure of directly applying the algorithm and the theoretical results of [34, Section 10.2].

Despite this, we can still design a two-layer nested optimization algorithm inspired by [34, Section 10.2]. Furthermore, our proposed algorithm can be proved to converge to the minimum CRLB under some necessary requirements, as will be explained in Section V-D.

The proposed algorithm is based on iterative maximization in the inner layer and the outer layer of (31). In the inner layer of (31), we use the stochastic Newton’s method [35] to update the estimate, given by

[TABLE]

where $b_{S,k}$ is the tracking step-size in Quasi-static Case, and $\boldsymbol{\varsigma}_{k}$ is the updating direction vector. This updating direction vector is a function of the observation vector $\textbf{y}_{k}$ and the latest estimated value of the channel parameter vector $\hat{\boldsymbol{\psi}}_{k-1}$, and is defined as below:

[TABLE]

And it is derived that

[TABLE]

where ${\textbf{e}}_{k}=\textbf{W}_{k}^{\text{H}}\textbf{a}\left(\hat{\textbf{x}}_{k-1}\right)$, $\hat{\textbf{y}}_{k}=\lvert\textbf{s}\rvert\hat{\beta}_{k-1}\textbf{W}_{k}^{\text{H}}\textbf{a}\left(\hat{\textbf{x}}_{k-1}\right)$, $\tilde{\textbf{e}}_{k1}=\hat{\beta}_{k-1}\textbf{W}_{k}^{\text{H}}\frac{\partial\textbf{a}\left(\hat{\textbf{x}}_{k-1}\right)}{\partial x_{1}}$, $\tilde{\textbf{e}}_{k2}=\hat{\beta}_{k-1}\textbf{W}_{k}^{\text{H}}\frac{\partial\textbf{a}\left(\hat{\textbf{x}}_{k-1}\right)}{\partial x_{2}}$ and $\hat{\textbf{V}}_{k}$ is given by

[TABLE]

In the outer layer of (31), assuming that the estimate of the channel parameter vector is accurate, i.e., ${{\boldsymbol{\psi}}=\hat{\boldsymbol{\psi}}_{k-1}}$, the EBM $\textbf{W}_{k}=\left[\textbf{w}_{k,1},\textbf{w}_{k,2},\textbf{w}_{k,3}\right]$ is obtained with

[TABLE]

Finally, the proposed tracking algorithm is summarized in Algorithm 1.

V-D Asymptotic Optimality Analysis

In this subsection, the convergence and the optimality of our proposed algorithm will be discussed. Since the entire proofs are very long, we will provide the main statements and ideas or clues here, and leave the proofs in the appendices.

The convergence and optimality will be stated in three steps as follows:

i) We prove that the estimate of the proposed tracking algorithm converges to a unique point with probability one given appropriate sequences of step-sizes.

ii) We prove that if the initial estimate is within the main lobe, i.e., $\hat{\textbf{x}}_{0}\in\mathcal{B}\left(\textbf{x}\right)$ and the step-size is appropriate, then the convergence point will be exactly the real channel parameter vector $\boldsymbol{\psi}$, with probability approaching one.

iii) Finally, if $\hat{\boldsymbol{\psi}}_{k}\to{\boldsymbol{\psi}}$ and the step-size is appropriate, then the tracking error of our algorithm converges to the minimum CRLB.

1) Convergence to a unique point

Since the observation vector $\textbf{y}_{k}$ is corrupted by the Gaussian noise vector, the updating direction vector $\boldsymbol{\varsigma}_{k}$ in (33) is also a random vector, and can be expressed as follows:

[TABLE]

where $\textbf{f}_{\boldsymbol{\psi}}\left(\hat{\boldsymbol{\psi}}_{k-1}\right)$ is the deterministic part of $\boldsymbol{\varsigma}_{k}$ defined as below:

[TABLE]

which is a function of $\hat{\boldsymbol{\psi}}_{k-1}$ that takes $\boldsymbol{\psi}$ as a parameter vector. The zero-mean random part of $\boldsymbol{\varsigma}_{k}$, i.e., $\hat{\textbf{z}}_{k}$, is given by

[TABLE]

The randomness of $\boldsymbol{\varsigma}_{k}$ might cause the proposed algorithm to diverge. However, if we adopt the diminishing step-size as that in [34, 36, 37], i.e.,

[TABLE]

where $K_{S,0}\geq 0$ and $\epsilon_{S}>0$, then some convergence property can be obtained, as described in the following theorem:

Theorem 1** (Convergence to a Unique Stable Point).**

If we adopt the iterative method in (37), (38) and $b_{S,k}$ is given by (42) with $\epsilon_{S}>0$ and $K_{S,0}\geq 0$, then $\hat{\boldsymbol{\psi}}_{k}$ converges to a unique stable point of $\textbf{f}_{\boldsymbol{\psi}}\left(\hat{\boldsymbol{\psi}}_{k-1}\right)$ with probability one.

A point $\hat{\boldsymbol{\psi}}_{k-1}$ is called a stable point of $\textbf{f}_{\boldsymbol{\psi}}\left(\hat{\boldsymbol{\psi}}_{k-1}\right)$ when it satisfies two conditions: 1) $\textbf{f}_{\boldsymbol{\psi}}\left(\hat{\boldsymbol{\psi}}_{k-1}\right)=\textbf{0}$ and 2) $\frac{\partial\textbf{f}_{\boldsymbol{\psi}}\left(\hat{\boldsymbol{\psi}}_{k-1}\right)}{\partial\hat{\boldsymbol{\psi}}_{k-1}^{\text{T}}}$ is negative definite. Hence, the stable points set is defined as below:

[TABLE]

where $\textbf{A}\prec\textbf{0}$ denotes that the matrix A is negative definite. In our problem, $\textbf{f}_{\boldsymbol{\psi}}\left(\hat{\boldsymbol{\psi}}_{k-1}\right)$ defined in (40) is given by

[TABLE]

Proof of Theorem 1. See Appendix D.$\hfill\blacksquare$

By Theorem 1, for the general step-size in (42), $\hat{\boldsymbol{\psi}}_{k}$ converges to a unique stable point in $\mathcal{S}$.

2) Convergence to the channel parameter vector $\boldsymbol{\psi}$

According to (43) and (44), it is easy to verify that the channel parameter vector ${\boldsymbol{\psi}}$ is a stable point by the following two points:

1. $\beta\textbf{W}_{k}^{\text{H}}\textbf{a}\left(\textbf{x}\right)=\hat{\beta}_{k-1}\textbf{e}_{k}$ in (44) when ${\hat{\boldsymbol{\psi}}_{k-1}={\boldsymbol{\psi}}}$. Hence, $\textbf{f}_{\boldsymbol{\psi}}\left({\boldsymbol{\psi}}\right)=\textbf{0}$;

2)\frac{\partial\textbf{f}_{\boldsymbol{\psi}}\left(\hat{\boldsymbol{\psi}}_{k-1}\right)}{\partial\hat{\boldsymbol{\psi}}_{k-1}^{\text{T}}}\big{|}_{\hat{\boldsymbol{\psi}}_{k-1}={\boldsymbol{\psi}}}=-\textbf{J}_{4} by derivation, where $\textbf{J}_{4}$ is the $4$-order identity matrix. Thus, \frac{\partial\textbf{f}_{\boldsymbol{\psi}}\left(\hat{\boldsymbol{\psi}}_{k-1}\right)}{\partial\hat{\boldsymbol{\psi}}_{k-1}^{\text{T}}}\big{|}_{\hat{\boldsymbol{\psi}}_{k-1}={\boldsymbol{\psi}}} is negative definite.

Therefore, ${\boldsymbol{\psi}}$ is a stable point, i.e., ${\boldsymbol{\psi}}\in\mathcal{S}$.

Other stable points in $\mathcal{S}$ correspond to the local optimal points of the beam and channel tracking problem, which are out of the main lobe $\mathcal{B}(\textbf{x})$ in (10). Except for the channel parameter vector ${\boldsymbol{\psi}}$, the antenna array gain of other stable points in $\mathcal{S}$ is quite low, resulting in low tracking accuracy. Unfortunately, the estimate of the DPV x may jump out of the main lobe in the tracking process and converge to other local optimal points due to the existence of observation noise. Hence, one key challenge is to ensure that the tracking algorithm converges to ${\boldsymbol{\psi}}$ rather than other stable points. Then we develop the following theorem to deal with this challenge:

Theorem 2** (Convergence to the DPV x).**

If we adopt the iterative method in (37), (38) and (i) the initial estimate of x is within the main lobe, i.e., $\hat{\textbf{x}}_{0}\in\mathcal{B}\left(\textbf{x}\right)$; (ii) $b_{S,k}$ is given by (42) with $\epsilon_{S}>0$, then there exist some $K_{S,0}\geq 0$ and $R>0$ such that

[TABLE]

Proof.

See Appendix E. ∎

We have assumed that the beam estimator in Fig. 1 can output an initial estimate $\hat{\textbf{x}}_{0}$ within the main lobe $\mathcal{B}\left(\textbf{x}\right)$. Under the condition $\hat{\textbf{x}}_{0}\in\mathcal{B}\left(\textbf{x}\right)$, Theorem 2 tells us the probability of $\hat{\textbf{x}}_{k}\to\textbf{x}$ is related to $\frac{\lvert\textbf{s}\rvert^{2}}{\epsilon_{S}^{2}\sigma_{z}^{2}}$. Hence, we can reduce the step-size or increase the transmit SNR $\frac{\lvert\textbf{s}\rvert^{2}}{\sigma_{z}^{2}}$ to make sure that $\hat{\textbf{x}}_{k}\to\textbf{x}$ approaching probability one.

According to Theorem 1, $\hat{\boldsymbol{\psi}}_{k}$ converges to a unique stable point corresponding to a local optimal point. Hence, this unique stable point will be exactly ${\boldsymbol{\psi}}$ when $\hat{\textbf{x}}_{k}\to\textbf{x}$, i.e., $\hat{\boldsymbol{\psi}}_{k}\to{\boldsymbol{\psi}}$.

3) Convergence with the minimum CRLB

Finally, the following theorem is developed to tell us the tracking error of the proposed algorithm:

Theorem 3** (Convergence to $\boldsymbol{\psi}$ with the minimum CRLB).**

If we adopt the iterative method in (37), (38) and (i) $\hat{\boldsymbol{\psi}}_{k}\to{\boldsymbol{\psi}}$; (ii) $b_{S,k}$ is given by (42) with $\epsilon_{S}=1$ and any $K_{S,0}\geq 0$, then $\hat{\textbf{h}}_{k}-\textbf{h}$ is asymptotically Gaussian and

[TABLE]

Proof.

See Appendix F. ∎

By Theorem 1, Theorem 2 and Theorem 3, if $\hat{\textbf{x}}_{0}\in\mathcal{B}\left(\textbf{x}\right)$ and we adopt the step-size $b_{S,k}$ in (42) with $\epsilon_{S}=1$ and $K_{S,0}\geq 0$, then the minimum CRLB is achieved asymptotically with high probability.

VI Recursive Beam Tracking for Dynamic Case I : Performance Bound, Convergence and Optimality

In Dynamic Case I, the channel gain changes fast while the beam direction changes slowly. We assume that the beam direction keeps static, i.e., $\textbf{x}_{k}=\textbf{x}=\left[x_{1},x_{2}\right]^{\text{T}}$. Hence, the antenna gain in the direction of the arriving path also keeps static, i.e., $\eta\left(\textbf{x}_{k}\right)=\eta\left(\textbf{x}\right)$. When the channel gain $\beta_{k}^{c}$ changes fast, it is very difficult to establish theorems of tracking the channel gain and beam direction simultaneously, as in Section V. Fortunately, acquiring the beam direction information is sufficient for alignment in mmWave mobile communication with analog beamforming. Hence, we only focus on beam direction tracking in Dynamic Case I.

The channel gains of adjacent ECCs in this section are assumed to be independent of each other. In addition, different distributions of the channel gain $\beta_{k}^{c}$ can lead to different suitable tracking strategies. The tracking strategy designed for one distribution of the channel gain may deteriorate sharply when applied to other distributions. Hence, each type of channel gain distribution deserves studying, of which Rayleigh fading channel is a special case that is easier to be analyzed. This special case happens when quite a number of rays existing in a cluster are indistinguishable. In this section, we choose Rayleigh fading channels to study for Dynamic Case I, i.e., \beta_{k}^{c}\sim\mathcal{CN}\big{(}0,\big{(}\sigma_{\beta}^{c}\big{)}^{2}\big{)}. Although the theoretical results in this section are only applicable for Rayleigh fading channels, the proposed algorithm is found robust for other types of time-varying channels according to the numerical results in Section IX-C.

When the channel gain $\beta_{k}^{c}$ is Gaussian distributed, the equivalent channel gain $\beta\left(\textbf{x}\right)=\eta\left(\textbf{x}\right)\beta_{k}^{c}$ also satisfies Gaussian distribution with the variance given below:

[TABLE]

Correspondingly, the observation vector $\textbf{y}_{k}$ in (9) satisfies Gaussian distribution for a given DPV x and EBM $\textbf{W}_{k}$, i.e., $\textbf{y}_{k}\sim\mathcal{CN}\left(\textbf{0},\boldsymbol{\Sigma}_{\textbf{y},k}\right)$, where $\boldsymbol{\Sigma}_{\textbf{y},k}$ is the covariance matrix of $\textbf{y}_{k}$ defined as follows:

[TABLE]

According to (48), we can obtain the determinant of $\boldsymbol{\Sigma}_{\textbf{y},k}$:

[TABLE]

Then the conditional probability density function of $\textbf{y}_{k}$ is given by

[TABLE]

The following structure of this section is similar to Section V: we first formulate the beam tracking problem and provide the lower bound of it. Then we develop a tracking algorithm and prove this algorithm can converge to the minimum CRLB.

VI-A Problem Formulation

Since we only track the beam direction in Dynamic Case I, the estimation function in (12) is reformulated as follows:

[TABLE]

Let $\Xi_{DI,k}=\left\{\textbf{F}_{k}^{c},\textbf{F}_{DI,k}^{e}\right\}$ denote a beam tracking scheme set in $k$-th ECC: based on historical observation vectors $\textbf{y}_{1},\cdots,\textbf{y}_{k-1}$ along with the corresponding EBMs $\textbf{W}_{1},\cdots,\textbf{W}_{k-1}$, choose an appropriate EBM $\textbf{W}_{k}$, apply it to obtain $\textbf{y}_{k}$ and make an estimation of the DPV x in $k$-th ECC by using all EBMs and observations available. Hence, in $k$-th ECC, the tracking problem is formulated as:

[TABLE]

where the constraint (53) ensures that $\hat{\textbf{x}}_{k}$ is an unbiased estimate of the DPV x.

Before providing a specific tracking algorithm, we will first explore the performance bound of the problem in (52).

VI-B Cramér-Rao Lower Bound of Tracking Error

We now perform some theoretical analysis on the beam tracking problem. Based on the CRLB theory in [33], we introduce the following lemma to obtain the lower bound of the tracking error:

Lemma 4**.**

In Dynamic Case I, given $\textbf{W}_{1},\cdots,\textbf{W}_{k}$, the MSE of the DPV estimate in (52) is lower bounded as follows:

[TABLE]

*where the Fisher information matrix * $\textbf{I}_{DI}\left(\textbf{x},{{\bf{W}}_{l}}\right)$ is given by

[TABLE]

and the $p$-th row, $j$-th column $\left(p=1,2;j=1,2\right)$ of $\textbf{I}_{DI}\left(\textbf{x},{{\bf{W}}_{l}}\right)$ is derived by (56)

with $\textbf{g}_{l}$, $\tilde{{g}}_{l,p}$ and $\textbf{G}_{l,p}$ defined below:

[TABLE]

Proof.

See Appendix G. ∎

The CRLB in (54) is a function of the EBMs $\textbf{W}_{1},\ldots,\textbf{W}_{k}$. Similar to that in Quasi-static Case, we consider the normalized CRLB (by multiplying $k$):

[TABLE]

By optimizing only one EBM W, we can further get the minimum CRLB, given by

[TABLE]

Solving problem (59) yields the optimal EBM \textbf{W}_{\small DI}^{*}=\big{[}\textbf{w}_{\small DI,1}^{*},\textbf{w}_{\small DI,2}^{*},\textbf{w}_{\small DI,3}^{*}\big{]} with

[TABLE]

where $\left\{\boldsymbol{\Delta}_{\small DI,1}^{*},\boldsymbol{\Delta}_{\small DI,2}^{*},\boldsymbol{\Delta}_{\small DI,3}^{*}\right\}$ denotes the optimal set of exploration offsets in Dynamic Case I.

VI-C Asymptotically Optimal Set of Exploration Offsets

In general, the CRLB in (59) is a function of a set of system parameters including the equivalent channel gain parameter $\sigma_{\beta}^{2}$, the DPV x and the array size $M,\,N$. Hence, the optimal set of 2D exploration offsets should also be a function of these parameters. Since it is very hard to obtain the expression of this optimal set, we adopt numerical search to deal with this issue. However, as many parameters in (59) may affect the optimal result, numerical search has to be reconducted for different parameter sets, resulting in high complexity.

Fortunately, through our investigation, some useful properties of the minimum CRLB and the optimal set of exploration offsets are given to simplify the numerical search, as described in the following lemma:

Lemma 5**.**

In Dynamic Case I, the minimum CRLB ${C}_{DI}^{\min}(\boldsymbol{\psi})$ and the optimal set of exploration offsets $\left\{\boldsymbol{\Delta}_{\small DI,1}^{*},\boldsymbol{\Delta}_{\small DI,2}^{*},\boldsymbol{\Delta}_{\small DI,3}^{*}\right\}$ have the following three properties:

1) ${C}_{DI}^{\min}(\boldsymbol{\psi})$ and $\left\{\boldsymbol{\Delta}_{\small DI,1}^{*},\boldsymbol{\Delta}_{\small DI,2}^{*},\boldsymbol{\Delta}_{\small DI,3}^{*}\right\}$ are invariant to the DPV x;

2) $\frac{\lvert\textbf{s}\rvert^{2}\sigma_{\beta}^{2}}{\sigma_{z}^{2}}{C}_{DI}^{\min}(\boldsymbol{\psi})$ converges to constant values as $\frac{\lvert\textbf{s}\rvert^{2}\sigma_{\beta}^{2}}{\sigma_{z}^{2}}\to+\infty$;

3) ${C}_{DI}^{\min}(\boldsymbol{\psi})$ converges as $\emph{M},\,\emph{N}\to+\infty$ and there exists a fixed set of exploration offsets that are unrelated to the array size and $\frac{\lvert\textbf{s}\rvert^{2}\sigma_{\beta}^{2}}{\sigma_{z}^{2}}$, denoted as $\left\{\widetilde{\boldsymbol{\Delta}}_{DI,1}^{*},\,\widetilde{\boldsymbol{\Delta}}_{DI,2}^{*},\widetilde{\boldsymbol{\Delta}}_{DI,3}^{*}\right\}$, such that

[TABLE]

where $\widetilde{\textbf{W}}_{DI}^{*}=[\tilde{\textbf{w}}_{DI,1}^{*},\tilde{\textbf{w}}_{DI,2}^{*},\tilde{\textbf{w}}_{DI,3}^{*}]$ is obtained with

[TABLE]

Proof.

See Appendix H. ∎

Lemma 5 reveals that $\left\{\boldsymbol{\Delta}_{\small DI,1}^{*},\boldsymbol{\Delta}_{\small DI,2}^{*},\boldsymbol{\Delta}_{\small DI,3}^{*}\right\}$ is only related to the array size $M,\,N$ and $\frac{\lvert\textbf{s}\rvert^{2}\sigma_{\beta}^{2}}{\sigma_{z}^{2}}$. Hence, the numerical search times can be reduced to one for a particular array size $M,\,N$ and a particular $\frac{\lvert\textbf{s}\rvert^{2}\sigma_{\beta}^{2}}{\sigma_{z}^{2}}$. Numerically, we find later that even if $\left\{\boldsymbol{\Delta}_{\small DI,1}^{*},\boldsymbol{\Delta}_{\small DI,2}^{*},\boldsymbol{\Delta}_{\small DI,3}^{*}\right\}$ may change for different array sizes and $\frac{\lvert\textbf{s}\rvert^{2}\sigma_{\beta}^{2}}{\sigma_{z}^{2}}$, $\left\{\widetilde{\boldsymbol{\Delta}}_{\small DI,1}^{*},\,\widetilde{\boldsymbol{\Delta}}_{\small DI,2}^{*},\widetilde{\boldsymbol{\Delta}}_{\small DI,3}^{*}\right\}$ can be used to take the place of $\left\{\boldsymbol{\Delta}_{\small DI,1}^{*},\boldsymbol{\Delta}_{\small DI,2}^{*},\boldsymbol{\Delta}_{\small DI,3}^{*}\right\}$ as long as the antenna size $M,N$ and $\frac{\lvert\textbf{s}\rvert^{2}\sigma_{\beta}^{2}}{\sigma_{z}^{2}}$ are sufficiently large. Therefore, the numerical search times is reduced to one in the end. Similar to that in Quasi-static Case, $\left\{\widetilde{\boldsymbol{\Delta}}_{DI,1}^{*},\,\widetilde{\boldsymbol{\Delta}}_{DI,2}^{*},\widetilde{\boldsymbol{\Delta}}_{DI,3}^{*}\right\}$ is called the asymptotically optimal set of exploration offsets in Dynamic Case I in this paper.

By numerical search in the main lobe in (10), one asymptotically optimal set of exploration offsets $\left\{\widetilde{\boldsymbol{\Delta}}_{\small DI,1}^{*},\,\widetilde{\boldsymbol{\Delta}}_{\small DI,2}^{*},\widetilde{\boldsymbol{\Delta}}_{\small DI,3}^{*}\right\}$ can be obtained in TABLE III and Fig. 5. With this set of exploration offsets, a general way to generate the EBM $\widetilde{\textbf{W}}_{\small DI}^{*}$ is obtained by (61) to achieve the minimum CRLB.

By adopting $\left\{\widetilde{\boldsymbol{\Delta}}_{\small DI,1}^{*},\,\widetilde{\boldsymbol{\Delta}}_{\small DI,2}^{*},\widetilde{\boldsymbol{\Delta}}_{\small DI,3}^{*}\right\}$ to smaller size antenna arrays when $\frac{\lvert\textbf{s}\rvert^{2}\sigma_{\beta}^{2}}{\sigma_{z}^{2}}=0\,\text{dB}$, we compare the minimum CRLB and the CRLB achieved by $\left\{\widetilde{\boldsymbol{\Delta}}_{\small DI,1}^{*},\,\widetilde{\boldsymbol{\Delta}}_{\small DI,2}^{*},\widetilde{\boldsymbol{\Delta}}_{\small DI,3}^{*}\right\}$ in TABLE III. As illustrated in Fig. 6, when antenna number $M=N\geq 8$, we can approach the minimum CRLB with a relative error less than $0.1\%$ by using $\left\{\widetilde{\boldsymbol{\Delta}}_{\small DI,1}^{*},\,\widetilde{\boldsymbol{\Delta}}_{\small DI,2}^{*},\widetilde{\boldsymbol{\Delta}}_{\small DI,3}^{*}\right\}$.

By applying $\left\{\widetilde{\boldsymbol{\Delta}}_{\small DI,1}^{*},\,\widetilde{\boldsymbol{\Delta}}_{\small DI,2}^{*},\widetilde{\boldsymbol{\Delta}}_{\small DI,3}^{*}\right\}$ to different $\frac{\lvert\textbf{s}\rvert^{2}\sigma_{\beta}^{2}}{\sigma_{z}^{2}}$ when $M=N=8$, we compare the minimum CRLB and the CRLB achieved by $\left\{\widetilde{\boldsymbol{\Delta}}_{\small DI,1}^{*},\,\widetilde{\boldsymbol{\Delta}}_{\small DI,2}^{*},\widetilde{\boldsymbol{\Delta}}_{\small DI,3}^{*}\right\}$ in TABLE III. As illustrated in Fig. 7, when $\frac{\lvert\textbf{s}\rvert^{2}\sigma_{\beta}^{2}}{\sigma_{z}^{2}}\geq 0\,\text{dB}$, we can approach the minimum CRLB with a relative error less than $0.1\%$ by using $\left\{\widetilde{\boldsymbol{\Delta}}_{\small DI,1}^{*},\,\widetilde{\boldsymbol{\Delta}}_{\small DI,2}^{*},\widetilde{\boldsymbol{\Delta}}_{\small DI,3}^{*}\right\}$.

As a conclusion, it is practical to apply this asymptotically optimal set of exploration offsets to any antenna array with $M=N\geq 8$, any channel gain with $\frac{\lvert\textbf{s}\rvert^{2}\sigma_{\beta}^{2}}{\sigma_{z}^{2}}\geq 0\text{dB}$ and any direction.

VI-D Recursive Beam Tracking with Asymptotic Optimality Analysis

For the Rayleigh fading channels, it is crucial to acquire the variance of the equivalent channel gain in (47), while it is hindered by the unknown antenna gain $\eta\left(\textbf{x}\right)$. Fortunately, the estimate of $\eta\left(\textbf{x}\right)$ and $\sigma_{\beta}^{2}$ can be seen as approximately accurate for a given antenna element pattern as the estimate of the DPV $\hat{\textbf{x}}_{k}$ approach x. Hence, we assume a perfectly-known $\sigma_{\beta}^{2}$ here to design the algorithm in Dynamic Case I. The deterioration of the tracking performance caused by the estimation error of the antenna gain will be evaluated in Section IX-D.

The proposed tracker is motivated by the following maximum likelihood problem:

[TABLE]

Similar to that in Section V, we propose a two-layer nested optimization algorithm to find the solution of (62). Finally, the proposed tracking algorithm is given in Algorithm 2.

We now perform the asymptotic optimality analysis. According to [34, 36, 37], the diminishing step-size is adopted as follows:

[TABLE]

where $K_{DI,0}\geq 0$ and $\epsilon_{{\small DI}}>0$. Then we can prove that if the initial estimate $\hat{\textbf{x}}_{0}$ is within the main lobe and $\epsilon_{{\small DI}}=1$, the proposed algorithm can converge to x with the minimum CRLB with high probability, i.e.,

[TABLE]

The proof is similar to that in Section V and the details are omitted here since nothing new is provided in the proof.

VII Joint Beam and Channel Tracking for Dynamic Case II

In Dynamic Case II where both the channel gain $\beta_{k}^{c}$ and the DPV $\textbf{x}_{k}$ change fast, the observation vector $\textbf{y}_{k}$ satisfies normal distribution with $\textbf{y}_{k}\sim\mathcal{CN}\left(\lvert\textbf{s}\rvert\beta\left(\textbf{x}_{k}\right)\textbf{W}_{k}^{\text{H}}\textbf{a}(\textbf{x}_{k}),\sigma_{z}^{2}\textbf{J}_{3}\right)$ for a given channel parameter vector $\boldsymbol{\psi}_{k}$ and EBM $\textbf{W}_{k}$. Hence, the conditional probability density function of the observation vector $\textbf{y}_{k}$ is given by

[TABLE]

Establishing theorems of tracking, as in Section V and Section VI, is very difficult in Dynamic Case II. Even if the theoretical analysis is not conducted in this section, we still provide a tracking algorithm in this section.

Inspired by the asymptotically optimal tracking algorithm in Section V and Section VI, we design a similar joint beam and channel tracking algorithm in Algorithm 3.

Different from the step-size in Quasi-static Case and Dynamic Case I, we adopt constant step-size in Dynamic Case II as the diminishing step-size cannot track the fast-changing $\textbf{x}_{k}$ and $\beta_{k}$. The constant step-size $b_{\small DII,k}$ will be specified later.

VIII Computational Complexity

In this section, we evaluate the computational complexity of the proposed tracking algorithms in Quasi-static Case, Dynamic Case I and Dynamic Case II. We focus on the complex arithmetic operations in the tracking stage including complex multiplication and division, while complex addition and subtraction are omitted since they require much fewer operations. It seems that Algorithm 1, Algorithm 2 and Algorithm 3 require a huge number of complex arithmetic operations due to the Fisher information matrix inversion in each ECC. However, most of these calculation work can be finished off-line, by which the complex operations are greatly reduced. The following lemma is proposed to tell us the specific computational complexity:

Lemma 6**.**

If the number of offline complex arithmetic operations is ignored since it is much smaller than the online ones as the tracking process lasts, then

1) for Algorithm 1 in Quasi-static Case and Algorithm 3 in Dynamic Case II, 45 complex arithmetic operations are required in each ECC;

2) for Algorithm 2 in Dynamic Case I, 28 complex arithmetic operations are required in each ECC.

Proof.

See Appendix I. ∎

According to Lemma 6, our algorithms can efficiently work without high complexity.

IX Numerical Results

In this section, some numerical results will be provided to verify the performance of our proposed tracking algorithms for Quasi-static Case, Dynamic Case I and Dynamic Case II. Based on the model in Section II, the parameters are set as: $M\!=\!N\!=\!8$, the antenna spacing $d_{1}\!=\!d_{2}\!=\!\frac{\lambda}{2}$, and the transmit SNR is $\frac{\lvert\textbf{s}\rvert^{2}}{\sigma_{z}^{2}}=0\,\text{dB}$. The antenna element pattern is based on the 3GPP model [38]. The vertical cut and the horizontal cut of the radiation power pattern (normalized with 0 dB at the central direction) for each element are given below:

[TABLE]

where $\theta_{\text{3\,dB}}=\frac{13\pi}{36}$ is the 3 dB beamwidth in the vertical direction, $\phi_{\text{3\,dB}}=\frac{13\pi}{36}$ is the 3 dB beamwidth in the horizontal direction and $\eta_{\max}=30\,\text{dB}$ is the maximum attenuation.

The combined radiation power pattern (normalized) of each antenna element in [38] is given by

[TABLE]

As can be seen in Fig. 8, we define two direction regions: central direction region ($\eta_{\text{dB}}\left(\theta,\phi\right)\geq-6\,\text{dB}$) and edge direction region ($\eta_{\text{dB}}\left(\theta,\phi\right)\leq-18\,\text{dB}$). We will separately evaluate the tracking performance in these two direction regions afterwards.

IX-A Reference Algorithms

Reference algorithms include the compressed sensing algorithm in [14], the 3GPP New Radio (NR) tracking algorithm in [16], the extended Kalman filter (EKF) algorithm in [15] and the recursive beam and channel tracking (RBCT) algorithm in [10]. For the compressed sensing algorithm in [14], we randomly choose phase shifts from $\left\{\pm{1},\pm{j}\right\}$ and then use the sparse recovery algorithm to estimate the DPV, where a discrete Fourier transform (DFT) dictionary with a size of 1024 is utilized. As for the 3GPP NR tracking algorithm, the last estimated beam direction and its adjacent beam directions are probed, then the system determines whether to switch the estimate in current ECC according to the strength of the received signals. For the EKF algorithm, we extend the method in [15] to the 2D array and the three exploring directions used in each ECC form a regular triangle within half the main lobe of the DPV estimate. The original RBCT algorithm in [10] is designed for 1D system and cannot directly support 2D tracking. To compare with our algorithms in the 2D case, we use two symmetrical explorations to track each dimension of the 2D beam for the RBCT algorithm.

For the initial beam estimation stage in Fig. 1, an exhaustive beam sweeping is conducted. Then an initial estimate is obtained by using the orthogonal matching pursuit method in [39]. This ensures that the initial estimate of the DPV, i.e., $\hat{\textbf{x}}_{0}$, is within the main lobe in (10).

In the tracking stage, three explorations are conducted in each ECC for all the algorithms. For the RBCT algorithm in [10], we use a buffer to store the received observations and update the estimate when receiving four new observations.

It is worth pointing out that the compressed sensing algorithm in [14] does not require an initial estimate. To ensure the fairness of all the algorithms, we compare the tracking performance versus the total number of explorations (i.e., assuming the compressed sensing algorithm uses the same total number of explorations) used both in the initial beam estimation stage and the tracking stage.

IX-B Performance bound of the tracking algorithms

As revealed in Section V-A and Section VI-B, the CRLB is a function of the adopted EBMs. Since the EBMs of our algorithms and the reference algorithms are quite different, the corresponding achieved CRLBs are also different. In this subsection, we will compare the CRLBs achieved by the EBMs of these different algorithms to verify the superiority of our optimal EBMs.

The fixed AoA ($\theta$,$\,\phi$) in Quasi-static Case and Dynamic Case I is chosen evenly and randomly in $\theta\in\left[-\frac{\pi}{6},\frac{\pi}{6}\right],\phi\in\left[\frac{\pi}{3},\frac{2\pi}{3}\right]$. The corresponding antenna gain of each element varies from -5.2 dB to 0 dB via (72). For Quasi-static case, the fixed channel gain $\beta^{c}$ is modeled as Rician fading with a K-factor $\kappa$=15 dB, according to the channel model in [40]. For Dynamic Case I, the channel gains between different ECCs are independent of each other and modeled as Rayleigh fading with $\frac{\lvert\textbf{s}\rvert^{2}\left(\sigma_{\beta}^{c}\right)^{2}}{\sigma_{z}^{2}}=0\,\text{dB}$. All simulation results in this section are averaged over 1000 random system realizations.

Fig. 9 and Fig. 10 shows the CRLBs achieved by the EBMs of different algorithms. It can be observed that our optimal EBMs both in Quasi-static Case and Dynamic Case I can achieve lower CRLB compared with the EBMs used by the reference algorithms, which results from the fact that the EBMs of our algorithms are carefully optimized and proven to be optimal in theory. Since the CRLB illustrates the performance bound, the other four algorithms cannot perform better than our algorithms in potential tracking accuracy.

IX-C Results of Tracking Accuracy

In this subsection, we will evaluate the tracking performance of our algorithms in the central direction region. The AoA ($\theta$,$\,\phi$) as defined in Section II is chosen evenly and randomly in $\theta\in\left[-\frac{\pi}{6},\frac{\pi}{6}\right],\phi\in\left[\frac{\pi}{3},\frac{2\pi}{3}\right]$. The corresponding antenna gain of each element varies from -5.2 dB to 0 dB via (72).

1. Quasi-static Case

The channel gain $\beta^{c}$ is modeled as Rician fading with a K-factor $\kappa$=15 dB, according to the channel model in [40]. The step-size is set as $b_{S,k}=\frac{1}{k}$.

As can be observed in Fig. 11,

the tracking accuracy of the compressed sensing algorithm and the 3GPP NR tracking algorithm gradually keeps steady with the increasing number of explorations, as both of the two algorithms are grids-of-beam based approaches and the performance is restricted by the codebook resolution. The EKF algorithm in Fig. 11 also shows a similar feature, which results from the fact that the original EKF algorithm itself cannot efficiently track static parameters [41]. Although the tracking error of the RBCT algorithm gradually reduces as the number of exploration increases, it cannot achieve the corresponding CRLB in Fig. 9. This is caused by tracking the horizontal and the vertical directions separately in the RBCT algorithm, causing loss as against joint tracking. Compared with the four reference algorithms, our proposed JBCT algorithm can approach the minimum CRLB quickly and achieve much lower tracking error.

2. Dynamic Case I

The channel gains between different ECCs are independent of each other and Rayleigh fading channels are adopted in each ECC with $\frac{\lvert\textbf{s}\rvert^{2}\left(\sigma_{\beta}^{c}\right)^{2}}{\sigma_{z}^{2}}=0\,\text{dB}$. The step-size is set as $b_{DI,k}=\frac{1}{k}$. As the RBCT algorithm in [10] does not support fast-fading channel tracking, it is excluded from the reference algorithms in Dynamic Case I.

Fig. 12

indicates that the DPV MSE of our proposed RBT algorithm can converge to the minimum CRLB if the antenna gain can be perfectly known. Even with the estimated antenna gain for tracking, our algorithm can also converge to the minimum CRLB and achieve much lower tracking error than other algorithms. Hence, the estimation error of the antenna gain has little influence on the tracking performance of our algorithm when the AoA is in the central direction region.

Further, we evaluate the robustness of the proposed algorithm for other types of time-varying channels. The channel gains between different ECCs are still independent of each other while Rician fading channels are adopted in each ECC with a K-factor $\kappa$=15 dB and $\frac{\lvert\textbf{s}\rvert^{2}\left(\sigma_{\beta}^{c}\right)^{2}}{\sigma_{z}^{2}}=0\,\text{dB}$, where \big{(}\sigma_{\beta}^{c}\big{)}^{2} denotes the average energy gain of the Rician fading channel. The algorithm designed for the Rayleigh fading channel is adopted to track the Rician fading channel here. It can be observed in Fig. 13 that our algorithm can still converge and achieve much lower tracking error than existing algorithms. These results show that the proposed algorithm in Dynamic Case I is robust to different time-varying channels as long as the variance of the channel gain is known.

3. Dynamic Case II

In Dynamic Case II, the initial AoA ($\theta_{0}$,$\,\phi_{0}$) as defined in Section II is chosen evenly and randomly in $\theta_{0}\in\left[-\frac{\pi}{6},\frac{\pi}{6}\right],\phi_{0}\in\left[\frac{\pi}{3},\frac{2\pi}{3}\right]$. The AoA ($\theta_{k}$,$\,\phi_{k}$) is modeled as a random walk process with return, i.e., $\theta_{k+1}=\theta_{k}+\varpi_{k}^{\theta}\Delta\theta_{k}$, $\phi_{k+1}=\phi_{k}+\varpi_{k}^{\phi}\Delta\phi_{k}$, where $\Delta\theta_{k},\Delta\phi_{k}\sim\mathcal{N}(0,\delta_{{A}}^{2})$, and $\varpi_{k}^{\theta},\varpi_{k}^{\phi}\in\left\{-1,1\right\}$ denote the rotation direction. The rotation direction $\varpi_{k}^{\theta},\varpi_{k}^{\phi}$ are chosen such that $\theta_{k}$ varies in $\left[-\frac{\pi}{6},\frac{\pi}{6}\right]$ and $\phi_{k}$ varies in $\left[\frac{\pi}{3},\frac{2\pi}{3}\right]$. The channel gain is modeled as a first-order Gaussian-Markov process, i.e., $\beta_{k+1}^{c}=\rho\beta_{k}^{c}+\gamma_{k}$, where $\gamma_{k}\sim\mathcal{CN}(0,1-\rho^{2})$. We adopt $\rho=0.995$ in simulation. As for the step-size, numerical results show that when $b_{DII,k}=0.7$, the joint beam and channel tracking algorithm can track beams with higher velocity. Hence, the step-size is set as a constant $b_{DII,k}=0.7$.

Fig. 14

indicates the proposed JBCT algorithm in Dynamic Case II can achieve higher tracking accuracy than the other four algorithms. In addition, if we set a tolerance error in Fig. 14, then our algorithm can support higher angular velocities.

IX-D The impact of the antenna pattern

In this subsection, we will evaluate the impact of the antenna pattern by setting the AoA in the edge direction region in Fig. 8. Other parameters are the same as the setting in Section IX-C.

1. Quasi-static Case

The AoA ($\theta$, $\phi$) as defined in Section II is chosen evenly and randomly in $\theta\in\left[\frac{\pi}{2}-\frac{\pi}{60},\frac{\pi}{2}\right],\phi\in\left[\pi-\frac{\pi}{60},\pi\right]$. The corresponding antenna gain of each element is -30 dB via (72).

In Fig. 15, it can be seen that our JBCT algorithm still outperforms existing algorithms when the AoA is in the edge direction region. Nevertheless, compared with Fig. 11, the minimum CRLB cannot be achieved any more and the performance of all algorithms deteriorate if the transmit power keeps unchanged, i.e., $\frac{\lvert\textbf{s}\rvert^{2}}{\sigma_{z}^{2}}=0\,\text{dB}$. This can be explained by the decrease of the equivalent SNR when the AoA is in the edge direction region. If we compensate the gain loss in the edge direction region by increasing the transmit power by 30 dB, i.e., $\frac{\lvert\textbf{s}\rvert^{2}}{\sigma_{z}^{2}}=30\,\text{dB}$, then Fig. 15 demonstrates the performance can be greatly improved and the minimum CRLB can be achieved again. Furthermore, we compare the tracking performance in two cases: 1) $\frac{\lvert\textbf{s}\rvert^{2}}{\sigma_{z}^{2}}=30\,\text{dB}$ when $\theta\in\left[\frac{\pi}{2}-\frac{\pi}{60},\frac{\pi}{2}\right],\phi\in\left[\pi-\frac{\pi}{60},\pi\right]$, where the antenna gain of each element is -30 dB; 2) $\frac{\lvert\textbf{s}\rvert^{2}}{\sigma_{z}^{2}}=0\,\text{dB}$ when $\theta\in\left[-\frac{\pi}{120},\frac{\pi}{120}\right],\phi\in\left[\frac{\pi}{2}-\frac{\pi}{120},\frac{\pi}{2}+\frac{\pi}{120}\right]$, where the antenna gain of each element can be seen as 0 dB. It can be observed in Fig. 15 that the performance in these two cases is almost the same. This shows that the deterioration of the tracking performance in the edge direction region only results from the decrease of the equivalent SNR.

2. Dynamic Case I

The AoA ($\theta$,$\,\phi$) as defined in Section II is chosen evenly and randomly in $\theta\in\left[\frac{\pi}{3},\frac{\pi}{2}\right],\phi\in\left[\frac{5\pi}{6},\pi\right]$. The corresponding antenna gain of each element varies from $-30\,\text{dB}$ to -20.4 dB via (72). As can be observed in Fig. 16, our RBT algorithm still outperforms existing algorithms for Rayleigh fading channel when the AoA is in the edge direction region. Since the equivalent SNR decreases, the proposed algorithm cannot converge to the minimum CRLB as before. If we compensate the gain loss in the edge direction region by increasing the transmit power by 30 dB, i.e., $\frac{\lvert\textbf{s}\rvert^{2}}{\sigma_{z}^{2}}=30\,\text{dB}$, then our algorithm can still converge to the minimum CRLB and achieve the same performance as that in Fig. 12 with the perfectly-known antenna gain.

With the estimated antenna gain, our RBT algorithm cannot converge to the minimum CRLB as before even when $\frac{\lvert\textbf{s}\rvert^{2}}{\sigma_{z}^{2}}=\,30\text{dB}$. This is caused by the larger slope in the edge direction region compared with the central edge direction region. Hence, a small estimation error of the AoA can result in a large deviation of the estimated equivalent channel gain parameter $\sigma_{\beta}^{2}$ in (47), leading to the non-convergence when using estimated antenna gain in Fig. 16.

3. Dynamic Case II

The initial AoA ($\theta_{0}$,$\,\phi_{0}$) as defined in Section II is chosen evenly and randomly in $\theta_{0}\in\left[\frac{\pi}{3},\frac{\pi}{2}\right],\phi_{0}\in\left[\frac{5\pi}{6},\pi\right]$. The rotation direction $\varpi_{k}^{\theta},\varpi_{k}^{\phi}$ are chosen such that $\theta_{k}$ varies in $\left[\frac{\pi}{3},\frac{\pi}{2}\right]$ and $\phi_{k}$ varies in $\left[\frac{5\pi}{6},\pi\right]$. The corresponding antenna gain of each element varies from $-30\,\text{dB}$ to -20.4 dB via (72). As can be seen in Fig. 17, all the algorithms cannot efficiently track the channels when the AoA varies in the edge direction region, since the equivalent SNR decreases sharply. If we compensate the gain loss in the edge direction region by increasing the transmit power by 30 dB, i.e., $\frac{\lvert\textbf{s}\rvert^{2}}{\sigma_{z}^{2}}=30\,\text{dB}$, then our JBCT algorithm can still achieve lower tracking error and faster tracking speed.

IX-E Computational Complexity

We then evaluate the computational complexity of our proposed algorithms. As can be seen in TABLE IV, our algorithms require fewer complex operations than other algorithms except 3GPP NR. Compared with 3GPP NR, the proposed algorithms can achieve much more accurate tracking without greatly increasing the computational complexity.

X Conclusion

This paper focuses on fast accurate beam and channel tracking for 2D phased antenna arrays. We first give the minimum exploration overhead of joint 2D tracking in theory. Then three tracking algorithms are developed according to different practical time-varying channel models.

In Quasi-static Case, the optimal exploration offsets are derived which are proved to a) be unrelated to the channel gain and the beam direction, b) be determined only by the array size, and c) approach constants as the array size goes to infinity. Also, a joint beam direction and channel gain tracking algorithm is proposed and the tracking error is proved to converge to the minimum CRLB.

In Dynamic Case I, an algorithm for beam only tracking is proposed, and it is proved to converge and achieve the minimum CRLB on the beam direction.

In Dynamic Case II, a joint tracking algorithm of beam direction and channel gain is proposed with faster and more accurate performance presented by simulation results.

This work is the first step to beam and channel tracking with 2D phased antenna arrays. In future work, we will further study the following problems: i) establishing the corresponding theorems in Dynamic Case II; ii) jointly tracking multipath channels; iii) tracking at both the transmitter and the receiver.

Appendix A Proof of Lemma 1

If the EBVs are of the steering vector forms, i.e., $\textbf{w}_{k,i}=\frac{1}{\sqrt{MN}}\textbf{a}\left(\boldsymbol{\omega}_{k,i}\right)$, where $\boldsymbol{\omega}_{k,i}=\left[\omega_{k,i1},\omega_{k,i2}\right]^{\text{T}}$ denotes the $i$-th exploring direction vector in $k$-th ECC, then the noiseless complex observation equation for the $i$-th observation is given in (73),

where Step (a) follows the definition of $y_{a}(\boldsymbol{\Delta})$:

[TABLE]

with $\boldsymbol{\Delta}\triangleq\left[\delta_{1},\delta_{2}\right]^{\text{T}}$. In our real tracking problem, the exploring direction vector $\boldsymbol{\omega}_{k,i}$ should be ensured within the main lobe of $\textbf{x}_{k}$ in (10), i.e., $\lvert{\omega}_{k,i1}-{x}_{k,1}\rvert<1$ and $\lvert{\omega}_{k,i2}-{x}_{k,2}\rvert<1$. Hence, we have that $y_{a}\left(\boldsymbol{\omega}_{k,i}-\textbf{x}_{k}\right)>0$.

The complex observation equation in (73) contains two real equations, i.e., an amplitude equation and a phase angle equation. Therefore, $q$ amplitude equations and $q$ phase angle equations can be obtained after $q$ observations. If we set the first observation as a reference, then we can obtain $(q-1)$ relative amplitude equations and $(q-1)$ relative phase angle equations from the remaining $(q-1)$ observations. However, these phase angle equations are not independent, as to be explained below.

From (73), we can obtain the phase angle equation:

[TABLE]

Thus the relative phase angle equation of the $i$-th observation $y_{k,i}$ regarding the first observation $y_{k,1}\,(i\neq 1)$ can be obtained as below:

[TABLE]

where $\omega_{k,{11}}-{\omega_{k,i1}}$ and $\omega_{k,{12}}-{\omega_{k,i2}}$ are determined by the exploring direction vectors and unrelated to the channel parameter vector $\boldsymbol{\psi}_{k}$. From (75), we can know that once the exploring directions are determined, the relative phase angles are known constants unrelated to $\boldsymbol{\psi}_{k}$. In other words, the relative phase angle equations as revealed in (75) cannot provide any information for estimating $\boldsymbol{\psi}_{k}$.

Following the conclusion above, we analyze the minimum exploration overhead in the following two cases:

1. If we want to obtain the unique solution of $\boldsymbol{\psi}_{k}$ within one ECC, at least 4 independent real equations with respect to $\boldsymbol{\psi}_{k}$ are needed since $\boldsymbol{\psi}_{k}$ contains four independent real variables (i.e., the real part $\text{Re}\left\{\beta\left(\textbf{x}_{k}\right)\right\}$, the imaginary part $\text{Im}\left\{\beta\left(\textbf{x}_{k}\right)\right\}$ of the equivalent channel gain $\beta_{k}$ and the two direction parameters ${x}_{k,1},{x}_{k,2}$). After $q$ explorations in each ECC, we can obtain $q$ independent amplitude equations and only 1 independent phase angle equation, which is $q+1$ independent real equations with respect to $\boldsymbol{\psi}_{k}$ in total. Hence, at least 3 explorations are needed to obtain 4 independent real equations and estimate 4 independent real variables of $\boldsymbol{\psi}_{k}$.

2. If we only want to obtain the unique solution of $\textbf{x}_{k}$ within one ECC, at least 2 independent real equations with respect to $\textbf{x}_{k}$ are needed since $\textbf{x}_{k}$ contains two independent real variables (i.e., two direction parameters ${x}_{k,1},{x}_{k,2}$). It seems that fewer explorations are sufficient. However, we cannot obtain any absolute amplitude and phase information with respect to $\textbf{x}_{k}$ from one observation in (73) since $\beta\left(\textbf{x}_{k}\right)$ is unknown. In addition, the relative phase angles are constants unrelated to $\textbf{x}_{k}$. Thus, the phase angle equations are useless for estimating $\textbf{x}_{k}$. After $q$ explorations in each ECC, we can obtain $q-1$ independent relative amplitude equations with respect to $\textbf{x}_{k}$ in total. Hence, at least 3 explorations are needed to obtain 2 independent real equations and estimate 2 independent real variables of $\textbf{x}_{k}$.

Therefore, the proof is completed.

Appendix B Proof of Lemma 2

In problem (18), the constraint (14) ensures that $\hat{\textbf{h}}_{k}$ is an unbiased estimate of h. Consider each element of the channel vector h, i.e., ${h_{mn}}({{\boldsymbol{\psi}}})=\beta{e^{j2\pi\left({\frac{{m-1}}{M}{x_{1}}+\frac{{n-1}}{N}{x_{2}}}\right)}}$. Immediately we have $\mathbb{E}\left[{h_{mn}}(\hat{{\boldsymbol{\psi}}}_{k})\right]={h_{mn}}({{\boldsymbol{\psi}}})$ since $\mathbb{E}\left[\hat{\textbf{h}}_{k}\right]=\textbf{h}$. According to Section 3.8 of [33], if a function $f\left(\hat{\boldsymbol{\psi}}\right)$ is an unbiased estimate of $f\left(\boldsymbol{\psi}\right)$, i.e., $\mathbb{E}\left[f(\hat{{\boldsymbol{\psi}}})\right]=f({{\boldsymbol{\psi}}})$, then we can obtain that

[TABLE]

where $\operatorname{Var}[f(\hat{\boldsymbol{\psi}})]$ denotes the variance of $f(\hat{\boldsymbol{\psi}})$ and $\textbf{I}(\boldsymbol{\psi})$ is the corresponding Fisher information matrix.

Combining (18) and (76), we have

[TABLE]

where Step $(a)$ is obtained by substituting (76) into (B) and Step $(b)$ is due to the definition of V in (21).

As for the Fisher information matrix in (2), we can obtain $\frac{\partial\text{log}\,p_{S}\left(\textbf{y}_{l}|\boldsymbol{\psi},\textbf{W}_{l}\right)}{\partial{\beta^{\text{re}}}}$ as follows:

[TABLE]

Similarly, $\frac{\partial\text{log}\,p_{S}\left(\textbf{y}_{l}|\boldsymbol{\psi},\textbf{W}_{l}\right)}{\partial{\beta^{\text{im}}}}$, $\frac{\partial\text{log}\,p_{S}\left(\textbf{y}_{l}|\boldsymbol{\psi},\textbf{W}_{l}\right)}{\partial{x_{1}}}$, and $\frac{\partial\text{log}\,p_{S}\left(\textbf{y}_{l}|\boldsymbol{\psi},\textbf{W}_{l}\right)}{\partial{x_{2}}}$ are given as

[TABLE]

Hence, the gradient of $\text{log}\,p_{S}\left(\textbf{y}_{l}|\boldsymbol{\psi},\textbf{W}_{l}\right)$ is obtained as follows:

[TABLE]

With the help of (80), we can obtain that

[TABLE]

Substituting (80) and (81) into (2), the Fisher information matrix is given as follows:

[TABLE]

where in Step $(c)$ we have used the following property of $\text{Re}\left\{\cdot\right\}$:

[TABLE]

with $\textbf{u},\,\textbf{v}$ denoting column vectors and $\bar{\textbf{u}}$ denoting the conjugate of u. Step $(d)$ is due to the exchangeability of $\mathbb{E}\left[\cdot\right]$ and $\text{Re}\left\{\cdot\right\}$:

[TABLE]

Step $(e)$ is due to the i.i.d. circularly symmetric complex Gaussian property of each element of $\textbf{z}_{l}$, which means that $\mathbb{E}\left[\textbf{z}_{l}\textbf{z}_{l}^{\text{H}}\right]=\sigma_{z}^{2}\textbf{J}_{3}$, where $\textbf{J}_{3}$ is the 3-order identity matrix. Step $(f)$ in (LABEL:eq_stepb) results from the property of complex Gaussian noise:

[TABLE]

Therefore, the Fisher information matrix is derived in (B) and Lemma 2 is proved in the end.

Appendix C Proof of Lemma 3

Lemma 3 is proved in three steps:

Step 1: We prove that ${C}_{S}^{\min}(\boldsymbol{\psi})$ and $\left\{\boldsymbol{\Delta}_{S,1}^{*},\boldsymbol{\Delta}_{S,2}^{*},\boldsymbol{\Delta}_{S,3}^{*}\right\}$ are unrelated to the equivalent channel gain $\beta$.

The basic method is block matrix inversion. We first rewrite the Jacobian matrix V in (21) as follows:

[TABLE]

where $\textbf{V}_{1}$ and $\textbf{V}_{2}$ are given by

[TABLE]

It is clear that both $\textbf{V}_{1}$ and $\textbf{V}_{2}$ are unrelated to $\beta$. Besides, we can obtain the following properties of $\textbf{V}_{1}$:

[TABLE]

where $\bar{\textbf{V}}_{1}$ denotes the conjugate of ${\textbf{V}}_{1}$.

With the help of the Jacobian matrix V in (86), the Fisher information matrix in (2) can be divided into four $2\times 2$ matrices as follows:

[TABLE]

where $\bar{\beta}$ denotes the conjugate of $\beta$ and A, B, D are defined as:

[TABLE]

with $\textbf{J}_{2}$ denoting the 2-order identity matrix. By combining (88) and (98), we can obtain the properties of B:

[TABLE]

By using the block matrix inversion method, the inverse of the Fisher information matrix in (C) is given by

[TABLE]

where ${{\textbf{I}}_{i{p_{1}}}}$ and ${{\textbf{I}}_{i{p_{2}}}}\left(\beta\right)$ are defined in (101) and (102):

[TABLE]

The middle part of ${\bf{I}}_{i{p_{2}}}$, i.e., $(\lvert\beta\rvert^{2}{\textbf{D}}-\text{Re}\left\{\bar{\beta}\textbf{B}^{\text{H}}\right\}{\bf{A}}{{}^{-1}}$ $\text{Re}\left\{\beta\textbf{B}\right\})$, can be rewritten as follows:

[TABLE]

where Step $(a)$ results from the properties of B in (99) and the definition of A in (98), Step $(b)$ is due to the fact that A is a real matrix and Step $(c)$ is due to the definition of $\textbf{I}_{s}$:

[TABLE]

Therefore, we can rewrite $\textbf{I}_{ip_{2}}$ in (102) as follows:

[TABLE]

By combining (V-A) and (100), we can obtain that

[TABLE]

where Step $(d)$ is by substituting (86) and (101) into (106). Since both $\textbf{V}_{1}$ in (87) and A in (98) are unrelated to the equivalent channel gain $\beta$, the first part of (106), i.e., $\operatorname{Tr}\left\{\textbf{A}^{-1}\textbf{V}_{1}^{\text{H}}\textbf{V}_{1}\right\}$ are unrelated to $\beta$. By substituting (86) and (105), we can obtain the second part of (106), i.e., $\operatorname{Tr}\left\{{{{\bf{I}}_{i{p_{2}}}}}\left(\beta\right)\textbf{V}^{\text{H}}\textbf{V}\right\}$ in (C),

where Step $(e)$ and Step $(f)$ follow the properties of the B in (99) and the definition of A in (98). It is clear that $\operatorname{Tr}\left\{{{{\bf{I}}_{i{p_{2}}}}}\left(\beta\right)\textbf{V}^{\text{H}}\textbf{V}\right\}$ is also unrelated to $\beta$ because none of the matrix $\textbf{A},\textbf{B},\textbf{V}_{1},\textbf{V}_{2}$ is related to $\beta$. Hence, ${C}_{S}(\boldsymbol{\psi},\textbf{W})$ is unrelated to $\beta$.

Since ${C}_{S}(\boldsymbol{\psi},\textbf{W})$ is unrelated to $\beta$, the minimum CRLB ${C}_{S}^{\min}(\boldsymbol{\psi})$ in (27) and the optimal EBM $\textbf{W}_{S}^{*}$ are also unrelated to $\beta$. Hence, the optimal set of exploration offsets $\left\{\boldsymbol{\Delta}_{S,1}^{*},\boldsymbol{\Delta}_{S,2}^{*},\boldsymbol{\Delta}_{S,3}^{*}\right\}$ is unrelated to the equivalent channel gain $\beta$.

Step 2: We prove that ${C}_{S}^{\min}(\boldsymbol{\psi})$ and $\left\{\boldsymbol{\Delta}_{S,1}^{*},\boldsymbol{\Delta}_{S,2}^{*},\boldsymbol{\Delta}_{S,3}^{*}\right\}$ are unrelated to the DPV x.

Consider the CRLB in (V-A). we will first prove that the Fisher information matrix $\textbf{I}_{S}(\boldsymbol{\psi},\textbf{W})$ is unrelated to the DPV x. Next, we will prove that $\textbf{V}^{\text{H}}\textbf{V}$ is also unrelated to x. Then it is clear that the minimum CRLB and the optimal set of exploration offsets $\left\{\boldsymbol{\Delta}_{S,1}^{*},\boldsymbol{\Delta}_{S,2}^{*},\boldsymbol{\Delta}_{S,3}^{*}\right\}$ are unrelated to x.

The Fisher information matrix in (C) tells us that only $\textbf{W}^{\text{H}}\textbf{V}$ may be related to x, which is given by

[TABLE]

with $\textbf{W}^{\text{H}}\textbf{a}\left(\textbf{x}\right)$, $\textbf{W}^{\text{H}}\frac{\partial\textbf{a}\left(\textbf{x}\right)}{\partial x_{1}}$ and $\textbf{W}^{\text{H}}\frac{\partial\textbf{a}\left(\textbf{x}\right)}{\partial x_{2}}$ expanded as follows:

[TABLE]

Since the EBVs are of the steering vector forms, i.e., $\textbf{w}_{i}=\frac{1}{\sqrt{MN}}\textbf{a}\left(\textbf{x}+\boldsymbol{\Delta}_{i}\right)$, where $\boldsymbol{\Delta}_{i}=\left[\delta_{i1},\delta_{i2}\right]^{\text{T}}$ denotes the $i$-th exploration offset, the elements of $\textbf{W}^{\text{H}}\textbf{a}\left(\textbf{x}\right)$ and $\textbf{W}^{\text{H}}\frac{\partial\textbf{a}\left(\textbf{x}\right)}{\partial x_{1}}$ can be written in (110) and (111).

[TABLE]

As shown in (110) and (111), both $\textbf{w}_{i}^{\text{H}}\textbf{a}\left(\textbf{x}\right)$ and $\textbf{w}_{i}^{\text{H}}\frac{\partial\textbf{a}\left(\textbf{x}\right)}{\partial x_{1}}$ are unrelated to the DPV x. Similarly, $\textbf{w}_{i}^{\text{H}}\frac{\partial\textbf{a}\left(\textbf{x}\right)}{\partial x_{2}}$ is also unrelated to x. Therefore, $\textbf{W}^{\text{H}}\textbf{V}$ in (108) is unrelated to x. Hence, the whole Fisher information matrix in (C) is invariant to x.

As for $\textbf{V}^{\text{H}}\textbf{V}$, we write it in (C),

which shows that $\textbf{V}^{\text{H}}\textbf{V}$ is unrelated to x.

Now it is clear that the CRLB in (V-A), i.e., ${C}_{S}(\boldsymbol{\psi},\textbf{W})$, is unrelated to x because both the Fisher information matrix $\textbf{I}_{S}(\boldsymbol{\psi},\textbf{W})$ and $\textbf{V}^{\text{H}}\textbf{V}$ are unrelated to x. Therefore, the minimum CRLB in (27) and the optimal set of exploration offsets $\left\{\boldsymbol{\Delta}_{S,1}^{*},\boldsymbol{\Delta}_{S,2}^{*},\boldsymbol{\Delta}_{S,3}^{*}\right\}$ are invariant to the DPV x.

Step 3: We prove that ${C}_{S}^{\min}(\boldsymbol{\psi})$ converges as $\emph{M},\,\emph{N}\to+\infty$ and

[TABLE]

Let us go into the asymptotic features of (V-A). According to (110) and (111), when the antenna number M, $\emph{N}\to+\infty$, the limit of the $i$-th ($i=1,2,3$) element of $\textbf{W}^{\text{H}}\textbf{a}\left(\textbf{x}\right)$, $\textbf{W}^{\text{H}}\frac{\partial\textbf{a}\left(\textbf{x}\right)}{\partial x_{1}}$ and $\textbf{W}^{\text{H}}\frac{\partial\textbf{a}\left(\textbf{x}\right)}{\partial x_{2}}$ in (109) are given as follows:

[TABLE]

where $\operatorname{Sa}\left[t\right]\triangleq\frac{\sin{t}}{t}$. Hence, each element of $\textbf{W}^{\text{H}}\textbf{V}/\sqrt{MN}$ in (108) converges when $M,\,N\to+\infty$, which results in that ${\textbf{I}}_{S}(\boldsymbol{\psi},\textbf{W})/{MN}$ in (C) also converges. The limit is defined as follows:

[TABLE]

The limit of $\textbf{V}^{\text{H}}\textbf{V}$ in (C) is given by

[TABLE]

By combining (114) and (C), we obtain the limit of $C_{S}(\boldsymbol{\psi},\textbf{W})$ in (V-A) as $M,\,N\to+\infty$:

[TABLE]

which reveals that the CRLB in (27), i.e., $C_{S}(\boldsymbol{\psi},\textbf{W})$, converges. Hence, the minimum CRLB $C_{S}^{\min}(\boldsymbol{\psi})$ also converges.

Let

[TABLE]

Then we have

[TABLE]

where Step (g) results from (117). On the other hand, we have

[TABLE]

where Step (h) results from (27). Hence, we can obtain that

[TABLE]

Therefore, Lemma 3 gets proved.

Appendix D Proof of Theorem 1

According to (32) and (39), the tracking procedure can be rewritten as

[TABLE]

And it can be derived in (41) that

[TABLE]

Since $\hat{\textbf{z}}_{k}\triangleq\left[\hat{z}_{k,1},\hat{z}_{k,2},\hat{z}_{k,3}\right]$ in (122) is composed of three i.i.d. circularly symmetric complex Gaussian random variables, the expectation of $\hat{\textbf{z}}_{k}$ is $\mathbb{E}\left[\hat{\textbf{z}}_{k}\right]=\textbf{0}$ and the covariance matrix is given in (123),

where Step $(a)$ is the result of (80) and Step $(b)$ follows the definition of the Fisher information matrix in (2).

Assume $\{\mathcal{G}_{k}:k\geq 0\}$ is an increasing sequence of $\sigma$-fields of $\{\hat{\boldsymbol{\psi}}_{0},\hat{\boldsymbol{\psi}}_{1},\hat{\boldsymbol{\psi}}_{2},\ldots\}$, i.e., $\mathcal{G}_{k-1}\!\subset\!\mathcal{G}_{k}$, where $\mathcal{G}_{0}\!\overset{\Delta}{=}\!\sigma(\hat{\boldsymbol{\psi}}_{0})$ and $\mathcal{G}_{k}\!\overset{\Delta}{=}\!\sigma(\hat{\boldsymbol{\psi}}_{0},\hat{\mathbf{z}}_{1},\ldots,\hat{\mathbf{z}}_{k})$ for $k\geq 1$. Because the $\hat{\mathbf{z}}_{k}$’s are composed of i.i.d. circularly symmetric complex Gaussian random variables with zero mean, $\hat{\mathbf{z}}_{k}$ is independent of $\mathcal{G}_{k-1}$, and $\hat{\boldsymbol{\psi}}_{k-1}\!\in\!\mathcal{G}_{k-1}$. Hence, we have

[TABLE]

for $k\geq 1$ and $\boldsymbol{\varsigma}_{k}=\textbf{f}_{\boldsymbol{\psi}}\left(\hat{\boldsymbol{\psi}}_{k-1}\right)+\hat{\textbf{z}}_{k}$ is also independent of $\mathcal{G}_{k-1}$.

Theorem 5.2.1 in [37, Section 5.2.1] gives the conditions that ensure $\hat{\boldsymbol{\psi}}_{k}$ converges to a unique point with probability one when there are several stable points. Next, we will prove that if the step-size $b_{S,k}$ is given by (42) with any $\varepsilon_{S}>0$ and $K_{S,0}\geq 0$, then the joint beam and channel tracking algorithm in (38) satisfies the corresponding conditions below:

Step-size requirements:

[TABLE]

• 2)

It is necessary to prove that $\sup\nolimits_{k}\mathbb{E}\left[\left\|\textbf{f}_{\boldsymbol{\psi}}\left(\hat{\boldsymbol{\psi}}_{k-1}\right)+\hat{\mathbf{z}}_{k}\right\|_{2}^{2}\right]<+\infty$.

From (121) and (123), we have

[TABLE]

where Step $(c)$ is due to (123) and that $\hat{\mathbf{z}}_{k}$ is independent of $\textbf{f}_{\boldsymbol{\psi}}\left(\hat{\boldsymbol{\psi}}_{k-1}\right)$.

From (44), we have

[TABLE]

As the Fisher information matrix is invertible, we get

[TABLE]

Besides, $\mathbf{W}_{k}\!=\!\left[\mathbf{w}_{k,1},\mathbf{w}_{k,2},\mathbf{w}_{k,3}\right]$, ${\mathbf{e}}_{k}\!=\!\mathbf{W}_{k}^{\text{H}}\mathbf{a}(\hat{\textbf{x}}_{k\!-\!1})$, $\tilde{\textbf{e}}_{k1}=\hat{\beta}_{k-1}\textbf{W}_{k}^{\text{H}}\frac{\partial\textbf{a}\left(\hat{\textbf{x}}_{k-1}\right)}{\partial x_{1}}$, $\tilde{\textbf{e}}_{k2}=\hat{\beta}_{k-1}\textbf{W}_{k}^{\text{H}}\frac{\partial\textbf{a}\left(\hat{\textbf{x}}_{k-1}\right)}{\partial x_{2}}$, hence we have

[TABLE]

[TABLE]

and

[TABLE]

for $i=1,2,3$ and all possible $\textbf{w}_{k,i}$ and x, where $\left[\delta_{k,i1},\delta_{k,i2}\right]^{\text{T}}=\boldsymbol{\omega}_{k,i}-\textbf{x}$. Thus we can get

[TABLE]

Combining (128) and (132), we have

[TABLE]

According to (128), it is clear that $\operatorname{Tr}\left\{\mathbf{I}(\hat{\boldsymbol{\psi}}_{k-1},\mathbf{W}_{k})^{-1}\right\}$ $<+\infty$. Then, we can get that

[TABLE]

• 3)

The function $\textbf{f}_{\boldsymbol{\psi}}\left(\hat{\boldsymbol{\psi}}_{k-1}\right)$ should be continuous with respect to $\hat{\boldsymbol{\psi}}_{k-1}$.

By using (44), we know that each element of $\textbf{f}_{\boldsymbol{\psi}}\left(\hat{\boldsymbol{\psi}}_{k-1}\right)$ is continuous with respect to $\hat{\boldsymbol{\psi}}_{k-1}=\left[\hat{\beta}_{k-1}^{\text{re}},\hat{\beta}_{k-1}^{\text{im}},\hat{x}_{k-1,1},\hat{x}_{k-1,2}\right]^{\text{T}}$. Therefore, $\textbf{f}_{\boldsymbol{\psi}}\left(\hat{\boldsymbol{\psi}}_{k-1}\right)$ is continuous with respect to $\hat{\boldsymbol{\psi}}_{k-1}$.

• 4)

Let $\boldsymbol{\mu}_{k}=\mathbb{E}\left[\left.\textbf{f}_{\boldsymbol{\psi}}\left(\hat{\boldsymbol{\psi}}_{k-1}\right)+\hat{\mathbf{z}}_{k}\right|\mathcal{G}_{k-1}\right]-\mathbf{f}\left(\hat{\boldsymbol{\psi}}_{k-1},\boldsymbol{\psi}\right)$. We need to prove that $\sum_{k=1}^{+\infty}\left\|b_{S,k}\boldsymbol{\mu}_{k}\right\|_{2}<+\infty$ with probability one.

From (D), we get $\boldsymbol{\mu}_{k}=\mathbf{0}$ for all $k\geq 1$. So we have $\sum_{k=1}^{+\infty}\left\|b_{S,k}\boldsymbol{\mu}_{k}\right\|_{2}=0<+\infty$ with probability one.

By Theorem 5.2.1 in [37], $\hat{\boldsymbol{\psi}}_{k}$ converges to a unique stable point within the stable points set with probability one.

Appendix E Proof of Theorem 2

Theorem 2 is proven in three steps:

Step 1: Two continuous processes based on the discrete process $\hat{\boldsymbol{\psi}}_{k}=[\hat{\beta}^{\text{re}}_{k},\hat{\beta}^{\text{im}}_{k},\hat{x}_{k,1},\hat{x}_{k,2}]^{\text{\emph{T}}}$ are established, i.e., $\bar{\boldsymbol{\psi}}(t)\!\overset{\Delta}{=}\![\bar{\beta}^{\text{re}}(t),\bar{\beta}^{\text{im}}(t),\bar{x}_{1}(t),\bar{x}_{2}(t)]^{\text{\emph{T}}}$ and $\tilde{\boldsymbol{\psi}}^{k}(t)\!\overset{\Delta}{=}\![\tilde{\beta}^{\text{re},k}(t),\tilde{\beta}^{\text{im},k}(t),\tilde{x}_{1}^{k}(t),\tilde{x}_{2}^{k}(t)]^{\text{\emph{T}}}$.

The discrete time parameters are defined as: $t_{0}\overset{\Delta}{=}0$, $t_{k}\overset{\Delta}{=}\sum_{l=1}^{k}b_{S,l}$, $k\geq 1$. The first continuous process $\bar{\boldsymbol{\psi}}(t),t\geq 0$ is constructed as the linear interpolation of the sequence $\hat{\boldsymbol{\psi}}_{k},k\geq 0$, where $\bar{\boldsymbol{\psi}}(t_{k})=\hat{\boldsymbol{\psi}}_{k},k\geq 0$. Therefore, $\bar{\boldsymbol{\psi}}(t)$ is given by

[TABLE]

The second continuous process $\tilde{\boldsymbol{\psi}}^{k}(t)$ is the solution of the following ordinary differential equation (ODE):

[TABLE]

for $t\in[t_{k},\infty)$, where $\tilde{\boldsymbol{\psi}}^{k}(t_{k})=\bar{\boldsymbol{\psi}}(t_{k})=\hat{\boldsymbol{\psi}}_{k},k\geq 0$. Thus, $\tilde{\boldsymbol{\psi}}^{k}(t)$ can be given as

[TABLE]

Step 2: By using the two continuous processes $\bar{\boldsymbol{\psi}}(t)$ and $\tilde{\boldsymbol{\psi}}^{k}(t)$ constructed in Step 1, a sufficient condition for the convergence of the discrete process $\hat{\textbf{x}}_{k}$ is provided here.

We first construct a time-invariant set $\mathcal{I}$ that includes the DPV x within the main lobe, i.e., $\textbf{x}\in\mathcal{I}\subset\mathcal{B}(\textbf{x})$. Define $\tilde{\textbf{x}}^{0}(t)\triangleq\left[\tilde{x}_{1}^{0}(t),\tilde{x}_{2}^{0}(t)\right]^{\text{T}}$ and denote $\hat{\textbf{x}}_{\text{b}}=\tilde{\textbf{x}}^{0}(t_{\text{b}})$ as the beam direction of the process $\tilde{\boldsymbol{\psi}}^{0}(t)$ that is closest to the boundary of the main lobe, which is given by333 The boundary of the set $\mathcal{B}(\textbf{x})$ is denoted by $\partial\mathcal{B}(\textbf{x})$.

[TABLE]

Then we pick $\delta$ such that

[TABLE]

where $\left\|\textbf{u}\right\|_{-\infty}=\underset{l=1,2}{\min}{\left[\textbf{u}\right]_{l}}$ denotes the minimum element of u. Note that when $t\geq t_{b}$, the solution $\tilde{\boldsymbol{\psi}}^{0}(t)$ of the ODE (136) will approach the real equivalent channel gain $\beta$ and DPV x monotonically as time $t$ increases. Hence, we construct the invariant set $\mathcal{I}$ in (140).

An example of the invariant set $\mathcal{I}$ is shown in Fig. 18.

Then, a sufficient condition will be established in Lemma 7 that ensures $\hat{\textbf{x}}_{k}\!\in\!\mathcal{I}~{}\text{for}~{}k\!\geq\!0$, and hence from Corollary 2.5 in [36], we can obtain that $\hat{\textbf{x}}_{k}$ converges to x. Before giving Lemma 7, let us provide some definitions first:

• •

Pick $T>0$ such that the solution $\tilde{\boldsymbol{\psi}}^{0}(t),t\geq 0$ of the ODE (136) with \tilde{\boldsymbol{\psi}}^{0}(0)=\big{[}\hat{\beta}^{\text{re}}_{0},\hat{\beta}^{\text{im}}_{0},\hat{x}_{{0},1}, \hat{x}_{{0},2}\big{]}^{\text{T}} satisfies $\inf_{\textbf{v}\in\partial\mathcal{B}}\left|\textbf{v}\!-\!\tilde{\textbf{x}}^{0}(t)\right|\geq 2\delta$ for $t\geq T$. Since when $t\geq t_{b}$, $\tilde{\textbf{x}}^{0}(t)$ will approach the DPV x monotonically as time $t$ increases, one possible $T$ is given by

[TABLE]

where $[\cdot]_{i}$ denotes the $i$-th element of the vector.

• •

Let $T_{0}\overset{\Delta}{=}0$ and $T_{l+1}\overset{\Delta}{=}\min\left\{t_{i}:t_{i}\geq T_{l}+T,i\geq 0\right\}$ for $l\geq 0$. Then $T_{l+1}-T_{l}\in[T,T+b_{S,1}]$ and $T_{l}=t_{\tilde{k}(l)}$ for some $\tilde{k}(l)\uparrow+\infty$, where $\tilde{k}(0)=0$. Let $\tilde{\boldsymbol{\psi}}^{\tilde{k}(l)}(t)$ denote the solution of ODE (136) for $t\in I_{l}\overset{\Delta}{=}\left[T_{l},T_{l+1}\right]$ with $\tilde{\boldsymbol{\psi}}^{\tilde{k}(l)}(T_{l})=\bar{\boldsymbol{\psi}}(T_{l})$, $l\geq 0$.

Hence, we can obtain the following lemma:

Lemma 7**.**

If $\underset{t\in I_{l}}{\sup}\left\|\bar{\textbf{x}}(t)-\tilde{\textbf{x}}^{\tilde{k}(l)}(t)\right\|_{2}\leq\delta$ for all $l\geq 0$, then $\hat{\textbf{x}}_{k}\in\mathcal{I}~{}\text{for all}~{}k\geq 0$.

Proof.

If $\underset{t\in I_{l}}{\sup}\left\|\bar{\textbf{x}}(t)-\tilde{\textbf{x}}^{\tilde{k}(l)}(t)\right\|_{2}\leq\delta$ for all $l\geq 0$, then $\underset{t\in I_{l}}{\sup}\left|\bar{x}_{1}(t)-\tilde{{x}}_{1}^{\tilde{k}(l)}(t)\right|\leq\delta$ and $\underset{t\in I_{l}}{\sup}\left|\bar{x}_{2}(t)-\tilde{{x}}_{2}^{\tilde{k}(l)}(t)\right|\leq\delta$.

According to Lemma 1 in [42], $\hat{{x}}_{k,1}\in\mathcal{I}~{}\text{for all}~{}k\geq 0$ and $\hat{{x}}_{k,2}\in\mathcal{I}~{}\text{for all}~{}k\geq 0$. Hence, $\hat{\textbf{x}}_{k}\in\mathcal{I}~{}\text{for all}~{}k\geq 0$. ∎

*Step 3: We will derive the probability lower bound for the condition in Lemma 7, which is also a lower bound for $P\left(\left.\hat{\textbf{x}}_{k}\!\rightarrow\!\textbf{x}\right|\hat{\textbf{x}}_{0}\!\in\!\mathcal{B}\left(\textbf{x}\right)\right)$. *

We will derive the probability lower bound for the condition in Lemma 7, which results in the following lemma:

Lemma 8**.**

If (i) the initial point satisfies $\hat{\textbf{x}}_{0}\in\mathcal{B}(\textbf{x})$, (ii) $b_{S,k}$ is given by (42) with any $\epsilon_{S}>0$, then there exist $K_{S,0}\geq 0$ and $R>0$ such that

[TABLE]

Proof.

See Appendix J. ∎

Finally, by applying Lemma 8 and Corollary 2.5 in [36], we can obtain

[TABLE]

which completes the proof of Theorem 2.

Appendix F Proof of Theorem  3

If the step-size $b_{S,k}$ is given by (42) with any $\varepsilon_{S}>0$ and $K_{S,0}\geq 0$, the sufficient conditions are provided by Theorem 6.6.1 [34, Section 6.6] to prove the asymptotic normality of $\sqrt{k}\left(\hat{\textbf{x}}_{k}-\textbf{x}\right)$, i.e., $\sqrt{k}\left(\hat{\textbf{x}}_{k}-\textbf{x}\right)\overset{d}{\rightarrow}\mathcal{N}\left(0,\Sigma_{\textbf{x}}\right)$. With the condition that $\hat{\boldsymbol{\psi}}_{k}\rightarrow\boldsymbol{\psi}$, we can prove that the beam and channel tracking algorithm satisfies the conditions above and obtain the variance $\boldsymbol{\Sigma}_{\textbf{x}}$ as follows:

Equation (121) is supposed to satisfy: (i) there exists an increasing sequence of $\sigma$-fields $\{\mathcal{F}_{k}:k\geq 0\}$ such that $\mathcal{F}_{l}\!\subset\!\mathcal{F}_{k}$ for $l\!<\!k$, and (ii) the random noise $\hat{\mathbf{z}}_{k}$ is $\mathcal{F}_{k}$-measurable and independent of $\mathcal{F}_{k-1}$.

As is shown in Appendix D, there exists an increasing sequence of $\sigma$-fields $\{\mathcal{G}_{k}:k\geq 0\}$, where $\boldsymbol{\varsigma}_{k}$ is measurable with respect to $\mathcal{G}_{k-1}$ and independent of $\mathcal{G}_{k-1}$.

• 2)

$\hat{\textbf{x}}_{k}$ should converge to x almost surely as $k\rightarrow+\infty$.

We assume that $\hat{\boldsymbol{\psi}}_{k}\rightarrow\boldsymbol{\psi}$, hence $\hat{\textbf{x}}_{k}$ converges to x almost surely when $k\rightarrow+\infty$.

• 3)

The stable condition:

In (44), we rewrite $\mathbf{f}_{\boldsymbol{\psi}}\left(\hat{\boldsymbol{\psi}}_{k-1}\right)$ as follows:

[TABLE]

where $\mathbf{D}_{1}$ is given by

[TABLE]

Then the stable condition is obtained that:

[TABLE]

which leads to $\varepsilon_{S}>\frac{1}{2}$.

• 4)

The noise vector $\hat{\mathbf{z}}_{k}$ satisfies:

[TABLE]

and

[TABLE]

Let

[TABLE]

where Step $(a)$ is obtained from (123).

By Theorem 6.6.1 [34, Section 6.6], we have

[TABLE]

where

[TABLE]

Due to that $\lim_{k\rightarrow+\infty}\sqrt{{(k+K_{S,0})}/{k}}=1$, we have

[TABLE]

if $k\rightarrow+\infty$. Thus, we can get

[TABLE]

By adopting $\epsilon_{S}=1$ in (148), we can obtain

[TABLE]

Since $\hat{\boldsymbol{\psi}}_{k}\to\boldsymbol{\psi}$ as $k\to+\infty$, $\hat{\textbf{h}}_{k}-\textbf{h}$ is linear to $\hat{\boldsymbol{\psi}}_{k}-\boldsymbol{\psi}$. Hence, $\hat{\textbf{h}}_{k}-\textbf{h}$ is also asymptotically Gaussian.

Combining (76), (150) and (27), we can conclude that

[TABLE]

Appendix G Proof of Lemma 4

In problem (52), the constraint (53) ensures that $\hat{\textbf{x}}_{k}$ is an unbiased estimate of x. According to Section 3.7 of [33], if $\hat{\textbf{x}}$ is an unbiased estimate of x, then we can obtain that

[TABLE]

where $\operatorname{Cov}\left(\hat{\textbf{x}}\right)$ denotes the covariance matrix of $\hat{\textbf{x}}$, $\textbf{I}\left(\textbf{x}\right)$ is the corresponding Fisher information matrix and $\textbf{A}\succeq\textbf{0}$ means that the matrix A is nonnegative definite. From (152), we can get that

[TABLE]

which implies that the diagonal elements of the matrix on the left side of ’$\succeq$’ are nonnegative because all matrices are $2\times 2$ in (153). Therefore, we obtain that

[TABLE]

i.e.,

[TABLE]

which yields the result of (54).

Now we try to obtain the Fisher information matrix in (55). According to (48), the determinant and the inverse of the covariance matrix can be written as follows:

[TABLE]

Based on the definition in (57), the determinant and the inverse of the covariance matrix in (156) can be rewritten as

[TABLE]

In addition, with the help of (50), we can obtain that

[TABLE]

where $\frac{\lvert\boldsymbol{\Sigma}_{\textbf{y},k}\rvert}{\partial{x}_{p}}$ and $\frac{\partial\boldsymbol{\Sigma}_{\textbf{y},k}^{-1}}{\partial{x}_{p}}$ are given by (159) according to (157):

[TABLE]

By combining (48), (50), (55) and (158), we can obtain the $p$-th row, $j$-th column element of the Fisher information below:

[TABLE]

Then we substitute (48), (49), (57), (157), (159) into (G), which yields the result of (56).

Finally, the proof of Lemma 4 is completed.

Appendix H Proof of Lemma 5

The proof of property 1) in Lemma 5 is similar to the proof of that in Lemma 3. Hence, we focus on the proof of property 2) and property 3) in Lemma 5.

Consider the $p$-th row, $j$-th column element of the Fisher information matrix $\textbf{I}_{DI}\left(\textbf{x},{{\bf{W}}_{k}}\right)$ in (56). We can rewrite it in (H),

where Step $(a)$ is obtained by substituting (156) into (H). When $\frac{\lvert\textbf{s}\rvert^{2}\sigma_{\beta}^{2}}{\sigma_{z}^{2}}\to+\infty$, we can obtain the element of $\textbf{I}_{DI}\left(\textbf{x},{{\bf{W}}_{k}}\right)$ in (162),

which reveals that $\frac{\sigma_{z}^{2}}{\lvert\textbf{s}\rvert^{2}\sigma_{\beta}^{2}}\textbf{I}_{DI}\left(\textbf{x},{{\bf{W}}_{k}}\right)$ converges as $\frac{\lvert\textbf{s}\rvert^{2}\sigma_{\beta}^{2}}{\sigma_{z}^{2}}\to+\infty$. Then the property 2) of Lemma 5 is proved.

Let us see the property 3) in Lemma 5. Similar to Step 3 in Appendix C, we can obtain that ${C}_{DI}^{\min}(\boldsymbol{\psi})$ converge as $M,N\to+\infty$ and

[TABLE]

According to (57) and (113), $\textbf{g}_{k}$ is $\Theta\left(\sqrt{MN}\right)$ while $\tilde{{g}}_{k,p}$ and $\textbf{G}_{k,p}$ are $\Theta\left(MN\right)$. Hence, $\sigma_{z}^{2}\operatorname{Tr}\left\{\textbf{G}_{k,p}\textbf{G}_{k,j}\right\}$ can be omitted since it is $\Theta\left({\left(MN\right)}^{2}\right)$ while other parts are $\Theta\left(\left(MN\right)^{\frac{5}{2}}\right)$. Then the $p$-th row, $j$-th column element of the Fisher information matrix in (H) can be rewritten in (H),

which reveals that $\left\{\widetilde{\boldsymbol{\Delta}}_{\small DI,1}^{*},\,\widetilde{\boldsymbol{\Delta}}_{\small DI,2}^{*},\widetilde{\boldsymbol{\Delta}}_{\small DI,3}^{*}\right\}$ is unrelated to $\frac{\lvert\textbf{s}\rvert^{2}\sigma_{\beta}^{2}}{\sigma_{z}^{2}}$.

Finally, the proof is completed.

Appendix I Proof of Lemma 6

We first analyze the computational complexity of Algorithm 1, which is composed of three steps:

Step 1: We evaluate the computational arithmetic operations of the Fisher information matrix inversion.

The Fisher information matrix is obtained as follows:

[TABLE]

where $\textbf{V}_{k}$ is given by

[TABLE]

with Step (a) resulting from the definition of $\textbf{V}_{k}^{1}$ and $\textbf{V}_{k}^{2}$:

[TABLE]

By combining (37), (110) and (111), we can obtain that $\textbf{W}_{k}^{\text{H}}\textbf{V}_{k}^{1}$ and $\textbf{W}_{k}^{\text{H}}\textbf{V}_{k}^{2}$ are determined matrices that remain unchanged for different ECCs, given by

[TABLE]

where both $\textbf{U}_{1}$ and $\textbf{U}_{2}$ can be obtained by offline calculation. Hence, we can rewrite the Fisher information matrix in (165) as:

[TABLE]

where $\bar{\beta}_{k-1}$ denotes the conjugate of $\hat{\beta}_{k-1}$ and $\widetilde{\textbf{A}}$, $\widetilde{\textbf{B}}$, $\widetilde{\textbf{D}}$ are defined as:

[TABLE]

Note that $\widetilde{\textbf{A}}$ is a diagonal matrix with the same diagonal elements and the block matrices $\widetilde{\textbf{A}}$, $\widetilde{\textbf{B}}$, $\widetilde{\textbf{D}}$ can all be obtained by offline calculation.

Similar to the derivation in (100), the inverse of the Fisher information matrix in (173) can be calculated by using the block matrix inversion method, given by

[TABLE]

where ${{\widetilde{\textbf{I}}}_{i{p_{1}}}}$ and ${{\widetilde{\textbf{I}}}_{i{p_{2}}}}\left(\hat{\beta}_{k-1}\right)$ are defined in (182) and (183)

[TABLE]

with $\widetilde{\textbf{I}}_{s}$ defined as follows:

[TABLE]

Since ${\widetilde{\textbf{A}}}{{}^{-1}}$ in (182) can be obtained by offline calculation, ${{\widetilde{\textbf{{I}}}}_{i{p_{1}}}}$ requires none online complex arithmetic operations.

As for ${{\widetilde{\textbf{{I}}}}_{i{p_{2}}}}$, we can further rewrite it as a block matrix:

[TABLE]

where the four block matrices are given by (190).

Since $\widetilde{\textbf{A}}$, $\widetilde{\textbf{B}}$, $\widetilde{\textbf{D}},\widetilde{\textbf{I}}_{s}$ can all be obtained by offline calculation, ${{\widetilde{\bf{I}}}_{i{p_{2}}}}^{12}\left(\hat{\beta}_{k-1}\right)$ and ${{\widetilde{\bf{I}}}_{i{p_{2}}}}^{22}\left(\hat{\beta}_{k-1}\right)$ only require 6 online complex arithmetic operations. In addition, ${{\widetilde{\bf{I}}}_{i{p_{2}}}}^{21}\left(\hat{\beta}_{k-1}\right)$ requires none online complex arithmetic operations as it can be obtained directly from ${{\widetilde{\bf{I}}}_{i{p_{2}}}}^{12}\left(\hat{\beta}_{k-1}\right)$. As for ${{\widetilde{\bf{I}}}_{i{p_{2}}}}^{11}\left(\hat{\beta}_{k-1}\right)$, we can convert it to (191)

. Finally, ${{\widetilde{\bf{I}}}_{i{p_{2}}}}^{11}\left(\hat{\beta}_{k-1}\right)$ requires none online complex arithmetic operations. Hence, the calculation of ${{\widetilde{\bf{I}}}_{i{p_{2}}}}\left(\hat{\beta}_{k-1}\right)$ in (190) requires none online complex arithmetic operations for ${{\widetilde{\bf{I}}}_{i{p_{2}}}}^{11}\left(\hat{\beta}_{k-1}\right)$, 6 online complex arithmetic operations for ${{\widetilde{\bf{I}}}_{i{p_{2}}}}^{12}$, none online complex arithmetic operations for ${{\widetilde{\bf{I}}}_{i{p_{2}}}}^{21}\left(\hat{\beta}_{k-1}\right)$, and 5 online complex arithmetic operations for ${{\widetilde{\bf{I}}}_{i{p_{2}}}}^{11}\left(\hat{\beta}_{k-1}\right)$, which are 11 complex arithmetic operations in total.

In the end, the calculation of $\textbf{I}_{S}\left(\hat{\boldsymbol{\psi}}_{k-1},\textbf{W}_{k}\right)^{-1}$ in (181) requires 11 online complex arithmetic operations.

Step 2: We evaluate the computational arithmetic operations of \frac{\partial\text{log}\,p_{S}\left(\textbf{y}_{k}|\boldsymbol{\psi},\textbf{W}_{k}\right)}{\partial\boldsymbol{\psi}}\bigg{|}_{\boldsymbol{\psi}=\hat{\boldsymbol{\psi}}_{k-1}}.

We write the \frac{\partial\text{log}\,p_{S}\left(\textbf{y}_{k}|\boldsymbol{\psi},\textbf{W}_{k}\right)}{\partial\boldsymbol{\psi}}\bigg{|}_{\boldsymbol{\psi}=\hat{\boldsymbol{\psi}}_{k-1}} as follows:

[TABLE]

where ${\textbf{e}}_{k}=\textbf{W}_{k}^{\text{H}}\textbf{a}\left(\hat{\textbf{x}}_{k-1}\right)$, $\hat{\textbf{y}}_{k}=\lvert\textbf{s}\rvert\hat{\beta}_{k-1}\textbf{W}_{k}^{\text{H}}\textbf{a}\left(\hat{\textbf{x}}_{k-1}\right)$, $\tilde{\textbf{e}}_{k1}=\hat{\beta}_{k-1}\textbf{W}_{k}^{\text{H}}\frac{\partial\textbf{a}\left(\hat{\textbf{x}}_{k-1}\right)}{\partial x_{1}}$, $\tilde{\textbf{e}}_{k2}=\hat{\beta}_{k-1}\textbf{W}_{k}^{\text{H}}\frac{\partial\textbf{a}\left(\hat{\textbf{x}}_{k-1}\right)}{\partial x_{2}}$. Since $\textbf{W}_{k}^{\text{H}}\textbf{a}\left(\hat{\textbf{x}}_{k-1}\right),\,\textbf{W}_{k}^{\text{H}}\frac{\partial\textbf{a}\left(\hat{\textbf{x}}_{k-1}\right)}{\partial x_{1}},\,\textbf{W}_{k}^{\text{H}}\frac{\partial\textbf{a}\left(\hat{\textbf{x}}_{k-1}\right)}{\partial x_{2}}$ can all be obtained by offline calculation, $\hat{\textbf{y}}_{k}$ requires $3$ online complex arithmetic operations, ${\textbf{y}}_{k}-\hat{\textbf{y}}_{k}$ requires none complex arithmetic operations and $\tilde{\textbf{e}}_{k1},\,\tilde{\textbf{e}}_{k2}$ both require 3 complex arithmetic operations. Together with the inner-product calculation in (192), the final number of online complex arithmetic operations is $18$.

Step 3: We evaluate the total computational arithmetic operations.

Considering the multiplication of $\textbf{I}_{S}\left(\hat{\boldsymbol{\psi}}_{k-1},\textbf{W}_{k}\right)^{-1}$ and \frac{\partial\text{log}\,p_{S}\left(\textbf{y}_{k}|\boldsymbol{\psi},\textbf{W}_{k}\right)}{\partial\boldsymbol{\psi}}\bigg{|}_{\boldsymbol{\psi}=\hat{\boldsymbol{\psi}}_{k-1}} (16 complex arithmetic operations), and the updating direction vector plus the previous estimate (none complex arithmetic operation), the final number of online complex arithmetic operations is 45 in each ECC.

By using a similar method, the total number of complex computational arithmetic operations for Algorithm 2 (Algorithm 3) is $28$ ($45$) in each ECC.

Therefore, Lemma 6 gets proved.

Appendix J Proof of Lemma 8

The following lemmas are introduced to prove Lemma 8.

Lemma 9** (Lemma 3 [42]).**

Given $T$ by (141) and

[TABLE]

If there exists a constant $C>0$, which satisfies

[TABLE]

for all $k\geq 0$ and $1\leq l\leq k_{T}$, then

[TABLE]

where $L$ and $C_{\mathbf{f}}$ are defined in (200) and (201) separately.

Lemma 10** (Lemma 4 [43]).**

If $\{M_{i}:i=1,2,\ldots\}$ satisfies that: (i) $M_{i}$ is Gaussian distributed with zero mean, and (ii) $M_{i}$ is a martingale in $i$, then

[TABLE]

for any $\eta>0$.

Lemma 11** (Lemma 5 [43]).**

If given a constant $C>0$, then

[TABLE]

is increasing for all $0<v<C$.

Let $\boldsymbol{\xi}_{0}\overset{\Delta}{=}\mathbf{0}$ and $\boldsymbol{\xi}_{k}\overset{\Delta}{=}\sum_{l=1}^{k}b_{S,l}\mathbf{\hat{z}}_{l}$, $k\geq 1$, where $\mathbf{\hat{z}}_{l}$ is given in (122). With (135) and (137), we have for $t_{k+l},1\leq l\leq k_{T}$,

[TABLE]

and

[TABLE]

where $\underline{v}\overset{\Delta}{=}\max\left\{t_{k}:t_{k}\leq v,k\geq 0\right\}$ for $v\geq 0$.

To bound $\int_{t_{k}}^{t_{k+l}}\left[\mathbf{f}_{\boldsymbol{\psi}}\left(\tilde{\boldsymbol{\psi}}^{k}(v)\right)-\mathbf{f}_{\boldsymbol{\psi}}\left(\tilde{\boldsymbol{\psi}}^{k}(\underline{v})\right)\right]dv$ on the RHS of (J), we obtain the Lipschitz constant of function $\mathbf{f}_{\boldsymbol{\psi}}(\mathbf{v})$ considering the first variable $\mathbf{v}$, given by

[TABLE]

Similar to (127), for any $t\geq t_{k}$, we can obtain that there exists a constant $0<C_{\mathbf{f}}<+\infty$ such that

[TABLE]

Hence, we have

[TABLE]

where Step $(a)$ is due to (200), Step $(b)$ is due to the definition in (137), and Step $(c)$ is due to (201). Then, by subtracting $\tilde{\boldsymbol{\psi}}^{k}(t_{k+l})$ in (J) from $\bar{\boldsymbol{\psi}}(t_{k+l})$ in (198) and taking norms, the following inequality can be obtained from (200) and (J) for $k\geq 0,1\leq l\leq k_{T}$:

[TABLE]

Applying Lemma 9 to (J) and letting

[TABLE]

yields

[TABLE]

where $C_{e}\overset{\Delta}{=}e^{L(T+b_{S,1})}$, and $c(k)\overset{\Delta}{=}\sum_{i>k}b_{S,i}^{2}$. Letting $k=\tilde{k}(l)$ in (204), we have $k+k_{T}=\tilde{k}(l+1)$ due to the definition of $T_{l+1}=t_{\tilde{k}(l+1)}$ in Step 2 of Appendix E and

[TABLE]

Suppose that the step size $\{b_{S,k}:k>0\}$ satisfies

[TABLE]

for $l\geq 0$. Given $\underset{t\in I_{l}}{\sup}\left\|\bar{\textbf{x}}(t)-\tilde{\textbf{x}}^{\tilde{k}(l)}(t)\right\|\!>\!\delta$, we can obtain from (205) and (206) that

[TABLE]

Then, we get

[TABLE]

where Step $(d)$ is due to the independence of noise, i.e., $\boldsymbol{\xi}_{p}-\boldsymbol{\xi}_{\tilde{k}(l)},\tilde{k}(l)\leq p\leq\tilde{k}(l+1)$ are independent of $\hat{\textbf{x}}_{k},0\leq k\leq\tilde{k}(l)$.

The lower bound of the probability that the sequence $\{\hat{\textbf{x}}_{k}:k\geq 0\}$ remains in the invariant set $\mathcal{I}$ is given by

[TABLE]

where Step $(e)$ is due to Lemma 7, Step $(f)$ is due to Lemma 4.2 in [36], and Step $(g)$ is due to (207). Let $\left\|\cdot\right\|_{\infty}$ denote the max-norm, i.e., $\left\|\mathbf{u}\right\|_{\infty}=\max_{l}|[\mathbf{u}]_{l}|$. Note that for $\mathbf{u}\in\mathbb{R}^{D}$, $\left\|\mathbf{u}\right\|_{2}\leq\sqrt{D}\left\|\mathbf{u}\right\|_{\infty}$. Hence we have

[TABLE]

With the increasing $\sigma$-fields $\{\!\mathcal{G}_{k}\!:\!k\!\geq\!0\!\}$ defined in Appendix D, we have for $k\geq 0$,

$\boldsymbol{\xi}_{k}\!=\!\sum_{l=1}^{k}b_{S,l}\hat{\mathbf{z}}_{l}\sim\mathcal{N}(0,\sum_{l=1}^{k}b_{S,k}^{2}\mathbf{I}_{S}(\hat{\boldsymbol{\psi}}_{l\!-\!1},\!\mathbf{W}_{l})^{-1})$,

• 2)

$\boldsymbol{\xi}_{k}$ is $\mathcal{G}_{k}$-measurable, i.e., $\mathbb{E}\left[\left.\boldsymbol{\xi}_{k}\right|\mathcal{G}_{k}\right]=\boldsymbol{\xi}_{k}$,

• 3)

$\mathbb{E}\left[\left\|\boldsymbol{\xi}_{k}\right\|^{2}_{2}\right]=\sum_{l=1}^{k}b_{S,k}^{2}\operatorname{tr}\left\{\mathbf{I}_{S}(\hat{\boldsymbol{\psi}}_{l\!-\!1},\!\mathbf{W}_{l})^{-1}\right\}<+\infty$,

• 4)

$\mathbb{E}\left[\left.\boldsymbol{\xi}_{k}\right|\mathcal{G}_{l}\right]=\boldsymbol{\xi}_{l}$ for all $0\leq l<k$.

Therefore, $\left[\boldsymbol{\xi}_{k}\right]_{j},j=1,2,3,4$ is a Gaussian martingale with respect to $\mathcal{G}_{k}$, and satisfies

[TABLE]

where C_{\mathbf{I}}\!\overset{\Delta}{=}\!\max_{s}\max_{i\geq 1}\frac{\lvert\textbf{s}\rvert^{2}}{\sigma_{z}^{2}}\big{[}\mathbf{I}(\hat{\boldsymbol{\psi}}_{i\!-\!1},\!\mathbf{W}_{i})^{-1}\big{]}_{j,j}. Let $\eta\!=\!\frac{\delta}{4C_{e}}$, M_{i}\!=\!\big{[}\boldsymbol{\xi}_{\tilde{k}(l)+i}\big{]}_{j}-\big{[}\boldsymbol{\xi}_{\tilde{k}(l)}\big{]}_{j},j\!=\!1,2,3,4 and $p={\tilde{k}(l+1)-\tilde{k}(l)}$ in Lemma 10, then from (J) and (J), we can obtain

[TABLE]

Combining (J), (J) and (J), we have

[TABLE]

To use Lemma 11, we assume that the step-size $b_{S,k}$ satisfies

[TABLE]

Then, from Lemma 11, we can obtain

[TABLE]

for $c(\tilde{k}(l))-c(\tilde{k}(l+1))<c(\tilde{k}(l))\leq c(0)$. Hence, we have

[TABLE]

As $C_{e}=e^{L(T+b_{S,1})}$, $c(0)=\sum_{i>0}b_{S,i}^{2}$, and $b_{S,k},T,L$ are given by (42), (141), (200) respectively, we can obtain

[TABLE]

In (J), $0<\delta<\inf_{\textbf{v}\in\partial\mathcal{B}}\left\|\textbf{v}-\hat{\textbf{x}}_{\text{b}}\right\|$, (206) and (213) should be satisfied, where a sufficiently large $K_{S,0}\geq 0$ can make both (206) and (213) true.

To ensure that $\hat{\textbf{x}}_{0}+b_{S,1}\left[\mathbf{f}_{\boldsymbol{\psi}}\left(\hat{\boldsymbol{\psi}}_{0}\right)\right]_{3,4}$ does not exceed the main lobe $\mathcal{B}(\textbf{x})$, i.e., the first step-size $b_{S,1}$ satisfies

[TABLE]

we can obtain the maximum $\epsilon_{S}$ as follows

[TABLE]

Hence, from (J), we have

[TABLE]

Combining (212), (214) and (217), yields

[TABLE]

which completes the proof.