Channel Estimation in mmWave Hybrid MIMO System via Off-Grid Dirichlet Kernels

TL;DR

This paper introduces a novel off-grid channel estimation method for mmWave hybrid MIMO systems using Dirichlet kernels, improving accuracy without increasing discretization complexity.

Contribution

It proposes greedy algorithms leveraging Dirichlet kernel structures in Fourier domain for off-grid parameter estimation in mmWave MIMO systems.

Findings

Proposed algorithms outperform standard OMP in off-grid scenarios.

Numerical results show reduced reconstruction errors with the new methods.

The approach effectively handles continuous angular parameters without dense discretization.

Abstract

In this paper, we tackle channel estimation in millimeter-wave hybrid multiple-input multiple-output systems by considering off-grid effects. In particular, we assume that spatial parameters can take any value in the angular domain, and need not fall on predefined discretized angles. Instead of increasing the number of discretized points to combat off-grid effects, we use implicit Dirichlet kernel structure in the Fourier domain, which conventional compressed sensing methods do not use. We propose greedy low-complexity algorithms based on orthogonal matching pursuit (OMP); our core idea is to traverse the Dirichlet kernel peak using estimates of the discrete Fourier transform. We demonstrate the efficacy of our proposed algorithms compared to standard OMP reconstruction. Numerical results show that our proposed algorithms obtain smaller reconstruction errors when off-grid effects are accounted for.

Equations

$$\mathbf{Y}={\mathbf{W}}^{\text{H}}\mathbf{HFs}+{\mathbf{N}}_W,$$

$$\mathbf{H}=\sum _{l=1}^L{\alpha }_l{\mathbf{a}}_{\text{UE}}({\theta }_l){\mathbf{a}}_{\text{BS}}^{\text{H}}({\varphi }_l),$$

$$\begin{array}{rl}& [{\mathbf{a}}_{\text{BS}}(\varphi ){]}_m=\frac{1}{\sqrt{M}}{e}^{j\frac{2\pi }{\lambda }{d}_{\text{BS}}(m-1)\sin (\varphi )},m\in [M],\\ & [{\mathbf{a}}_{\text{UE}}(\theta ){]}_n=\frac{1}{\sqrt{N}}{e}^{j\frac{2\pi }{\lambda }{d}_{\text{UE}}(n-1)\sin (\theta )},n\in [N],\text{vspace}-.25cm\end{array}$$

$$\begin{array}{rl}& \mathbf{\Phi }=\{\hat{\varphi }:\frac{1+\sin (\hat{\varphi })}{2}=\frac{m-1}{M};m\in [G_{\text{BS,D}}]\},\\ & \mathbf{\Theta }=\{\hat{\theta }:\frac{1+\sin (\hat{\theta })}{2}=\frac{n-1}{N};n\in [G_{\text{UE,D}}]\},\end{array}$$

$$\begin{array}{rl}& [{\mathbf{A}}_{\text{BS,D}}{]}_{m,m^{\prime }}=\frac{1}{\sqrt{M}}{e}^{j2\pi (m-1)(\frac{m^{\prime }-1}{M}-\frac{1}{2})};m,m^{\prime }\in [M],\\ & [{\mathbf{A}}_{\text{UE,D}}{]}_{n,n^{\prime }}=\frac{1}{\sqrt{N}}{e}^{j2\pi (n-1)(\frac{n^{\prime }-1}{N}-\frac{1}{2})};n,n^{\prime }\in [N].\end{array}$$

$$\begin{array}{rl}{\mathbf{H}}_{\text{V}}& ={\mathbf{A}}_{\text{UE,D}}^{\text{H}}\mathbf{H}{\mathbf{A}}_{\text{BS,D}}=\sum _{l=1}^L{\alpha }_l{\hat{\mathbf{a}}}_{\text{UE}}({\theta }_l){\hat{\mathbf{a}}}_{\text{BS}}^{\text{H}}({\varphi }_l),\end{array}$$

$$\displaystyle\begin{aligned} \big{[}\mathbf{\hat{a}}_{\text{UE}}(\theta_{l})\big{]}_{n^{\prime}}&=\frac{1}{{N}}\sum_{i=0}^{N-1}e^{j2\pi i(\frac{n^{\prime}-1}{N}-\frac{1}{2})}e^{j\pi i\sin{\theta_{l}}}\\ &=\frac{1}{{N}}\frac{\sin(\pi\vartheta_{n^{\prime},l}N)}{\sin{(\pi\vartheta_{n^{\prime},l})}}e^{-j\pi\vartheta_{n^{\prime},l}(N-1)},\end{aligned}$$

$$\displaystyle\begin{aligned} \big{[}\mathbf{H}_{\text{V}}\big{]}_{m^{\prime},n^{\prime}}=\sum_{l=1}^{L}\alpha_{l}\mathcal{D}(\varphi_{m^{\prime},l},\vartheta_{n^{\prime},l})\frac{e^{-j\pi\vartheta_{n^{\prime},l}(N-1)}}{e^{-j\pi\varphi_{m^{\prime},l}(M-1)}},\end{aligned}$$

$$\begin{array}{r}{}_{m^{\star },n^{\star }}=\sum _{l=1}^L{\alpha }_l\mathcal{D}({\phi }_{m^{\star },l},{\vartheta }_{n^{\star },l})\frac{e^{-j\pi {\vartheta }_{n^{\star },l}(N-1)}}{e^{-j\pi {\phi }_{m^{\star },l}(M-1)}}.\end{array}$$

$$\mathbf{y}=\mathbf{A}\text{vec}([{\mathbf{H}}_{\text{V}}{]}_{m^{\prime },n^{\prime }})+{\mathbf{n}}_W,$$

$$\begin{array}{rl}{P}_K& =\frac{4}{MN}\sum _{i=1}^{K/4}\sum _{j=1}^{K/4}{\left|\frac{\sin (M\pi \frac{(2i-1)}{2})}{\sin (\pi \frac{(2i-1)}{2})}\frac{\sin (N\pi \frac{(2j-1)}{2})}{\sin (\pi \frac{(2j-1)}{2})}\right|}^2,\\ P_T& =\frac{4}{MN}\sum _{i=1}^{M/4}\sum _{j=1}^{N/4}{\left|\frac{\sin (M\pi \frac{(2i-1)}{2})}{\sin (\pi \frac{(2i-1)}{2})}\frac{\sin (N\pi \frac{(2j-1)}{2})}{\sin (\pi \frac{(2j-1)}{2})}\right|}^2.\end{array}$$

$$\begin{array}{rl}& {\widehat{{\alpha }_l}}^{\star },{\hat{m}}_l^{\star },{\hat{n}}_l^{\star }=\underset{{\alpha }_l^{\star },m_l^{\star },n_l^{\star }}{\min }||\mathbf{y}-\mathbf{A}\text{vec}([{\mathbf{H}}_{\text{V}}{]}_{m^{\prime },n^{\prime }})|{|}_2^2\end{array}$$

$$\begin{array}{rl}& \text{s.t.}1\le m,m_l^{\star }\le M,1\le n,n_l^{\star }\le N,m^{\prime }\in [M],n^{\prime }\in [N],\\ & {\mathbf{H}}_{\text{V}}=\sum _{l=1}^L\frac{{\alpha }_l^{\star }}{MN}\frac{\sin (\pi (m-m_l^{\star }))\sin (\pi (n-n_l^{\star }))e^{-j\pi (m-m_l^{\star })}}{\sin (\frac{\pi }{M}(m-m_l^{\star }))\sin (\frac{\pi }{N}(n-n_l^{\star }))e^{-j\pi (n-n_l^{\star })}}.\end{array}$$

$$m_l^{\star }=m^{\prime }+\frac{\tan (\frac{\pi }{M})}{\frac{\pi }{M}}{\Delta }_m;n_l^{\star }=n^{\prime }+\frac{\tan (\frac{\pi }{N})}{\frac{\pi }{N}}{\Delta }_n$$

$${\Delta }_m=\mathbb{R}\left(\frac{[{\mathbf{H}}_{\text{V}}{]}_{m^{\prime }-1,n^{\prime }}-[{\mathbf{H}}_{\text{V}}{]}_{m^{\prime }+1,n^{\prime }}}{2[{\mathbf{H}}_{\text{V}}{]}_{m^{\prime },n^{\prime }}-[{\mathbf{H}}_{\text{V}}{]}_{m^{\prime }-1,n^{\prime }}-[{\mathbf{H}}_{\text{V}}{]}_{m^{\prime }+1,n^{\prime }}}\right)$$

$${\Delta }_n=\mathbb{R}\left(\frac{[{\mathbf{H}}_{\text{V}}{]}_{m^{\prime },n^{\prime }-1}-[{\mathbf{H}}_{\text{V}}{]}_{m^{\prime },n^{\prime }+1}}{2[{\mathbf{H}}_{\text{V}}{]}_{m^{\prime },n^{\prime }}-[{\mathbf{H}}_{\text{V}}{]}_{m^{\prime },n^{\prime }-1}-[{\mathbf{H}}_{\text{V}}{]}_{m^{\prime },n^{\prime }+1}}\right)$$

$${\alpha }_l^{\star },{\hat{m}}_l^{\star },{\hat{n}}_l^{\star }=\underset{{\alpha }_l,m_l^{\star },n_l^{\star }}{\min }||{\mathbf{y}}_{\text{res}}-\mathbf{A}\text{vec}([{\mathbf{H}}_{{\text{V}}_l}{]}_{m^{\prime },n^{\prime }})|{|}_2^2$$

$$\text{s.t.}m\in (m^{\prime }-1,m^{\prime }+1),n\in (n^{\prime }-1,n^{\prime }+1),$$

$$\text{HVl=MNαlsin(Mπ(m−ml⋆)sin(π(m−ml⋆))sin(Nπ(n−nl⋆))sin(π(n−nl⋆))e−jπ(m−ml⋆)e−jπ(n−nl⋆)}.$$

Full text

I Introduction and Related Work

One of the most promising features of next-generation wireless systems is to use high-frequency high-bandwidth signals in millimeter-wave (mmWave) frequency bands. These mmWave bands combined with multiple-input multiple-output (MIMO) technology have great potential in delivering higher data rates, higher spectral efficiency, and lower latency, exceeding the performance of traditional cellular systems operating at sub-6 GHz bands. Conventional mmWave MIMO architectures use a large number of antennas, which results in high cost and power consumption, making it difficult to assign a radio frequency (RF) chain per antenna. To curtail these issues, a hybrid analog/digital beamforming (HADB) architecture is adapted at mmWave bands [1, 2].

The HADB architecture complicates the channel estimation process, because only the low dimensional signals pre-combined by the analog combiner are available at baseband, which severely degrades the channel estimation process. The accuracy with which the channel is estimated plays a critical role in physical layer performance as it directly affects receiver design, e.g., channel equalization [3] and radio resource management [4]. To overcome these challenges, channel estimation algorithms based on compressed sensing (CS) [5, 2] have been proposed. These CS-based methods are based on virtual channel models [6], which provide a virtual angular representation of MIMO channels.

The virtual channel model describes the channel with respect to (w.r.t.) fixed basis functions corresponding to spatial angles within a finite discrete dictionary. In other words, the continuous parameter space of spatial angular features is discretized into a finite set of pre-defined spatial angles, which emphasizes the sparse representation of the MIMO channels. The estimation accuracy of CS methods based on this discretization is limited by the number of points in the dictionary. Although this discretization procedure yields state-of-art performance, it has several intrinsic disadvantages [7], including the off-grid effect.

A natural yet inefficient approach to reduce off-grid effects is to increase the number of discretized points, corresponding to increased grid resolution. This approach not only increases the mutual coherence of the dictionary matrix, leading to loss of the restricted isometric property, but also increases the problem dimension, which requires more computation [8]. An alternative is to tackle off-grid effects upfront without increasing the grid size. For example, in the context of channel estimation, Gurbuz et al. [8] provide a controlled perturbation mechanism for spatial angular parameters based on orthogonal matching pursuit (OMP) [9]. Other related works involve an improved off-grid sparse Bayesian algorithm [10], and a grid-less CS technique developed via atomic norm minimization [11]. Although these methods all tackle off-grid issues, they are computationally prohibitive, which motivates us to develop and analyze robust low-complexity channel estimation algorithms that account for off-grid effects.

Interestingly, standard CS methods based on sparsity fail to leverage Dirichlet structure in the Fourier domain. We exploit this structure to improve the channel estimation process. In particular, we propose low-complexity algorithms based on OMP [9], owing to its computational tractability. Our numerical results show that while accounting for off-grid effects, our proposed algorithms obtain smaller channel reconstruction errors compared to standard OMP algorithms.

Notation: Vectors and matrices are represented by lower-case and capital boldface letters, respectively (e.g.: a and A). The transpose, conjugate, conjugate transpose, and pseudo-inverse of a matrix A are denoted by $\mathbf{A}^{\text{T}}$, $\mathbf{A}^{\text{H}}$, $\mathbf{A}^{*}$, and $\mathbf{A}^{\dagger}$, respectively. For a non-negative integer $K$, we denote the set $\{1,2,\ldots,K\}$ by $[K]$; $\otimes$ denotes the Kronecker product; vec($\mathbf{A}$) denotes the vectorized version of the matrix $\mathbf{A}$; $\mathbb{R}$ is the real part of a complex number; and $\min(a,b)$ is the minimum of the scalars $a$ and $b$.

II System Model and Channel Model

II-A System Model

Consider a mmWave MIMO network comprised of a base station (BS) communicating with generic user equipment (UE), both equipped with a uniform linear array (ULA). We assume the BS is equipped with $M$ antennas, $M_{\text{RF}}$ RF chains, and $M_{\text{DS}}$ data streams. Similarly, the UE is assumed to be equipped with $N$ antennas, $N_{\text{RF}}$ RF chains, and $N_{\text{DS}}$ data streams. Typically, it is assumed that $M_{\text{DS}}\leq M_{\text{RF}}\leq M$ and $N_{\text{DS}}\leq N_{\text{RF}}\leq N$. With the HADB MIMO processing structure [1], the received signal at the UE is expressed,

[TABLE]

where $\mathbf{Y}$ $\in\mathcal{C}^{N_{\text{t}}\times M_{\text{t}}}$ is the received measurement matrix at the UE, and it is assumed that the UE uses $N_{\text{t}}$ combiners for each $M_{\text{t}}$ beamforming vector used by the BS. The received measurements at the UE for each beamforming vector are arranged in columns. The matrices $\mathbf{H}$ $\in\mathcal{C}^{N\times M}$, $\mathbf{W}$ $\in\mathcal{C}^{N\times N_{\text{t}}}$, and $\mathbf{F}$ $\in\mathcal{C}^{M\times M_{\text{t}}}$ represent the channel matrix from the BS to UE, the combined effect of the RF/baseband combiner, and precoder matrices, respectively. The noise matrix at the UE after the combiner operation is $\mathbf{N}_{W}=\mathbf{W}^{\text{H}}[\mathbf{n}_{1},\ldots,\mathbf{n}_{M_{\text{t}}}]$, where $\mathbf{n}_{\text{i}}\in\mathcal{C}^{N\times 1};\forall\text{i}\in[M_{\text{t}}]$, follows a circularly symmetric independent and identically distributed Gaussian distribution, $\mathcal{CN}(0,\sigma_{n}^{2})$ with noise variance $\sigma_{n}^{2}$. Further, $\mathbf{s}$ is assumed to be known at the BS and omitted hereafter.

II-B Channel Model

Based on [12], the mmWave channel model is given by:

[TABLE]

where $\alpha_{l}$ is the complex gain associated with $l^{\text{th}}$ multi-path component (MPC) between the BS and UE. The number of MPCs is $L$ with $L\,{\ll}\,{M;N}$ may itself be time-varying due to the mobility of the UE and the surrounding scatterers [13]. The terms $\mathbf{a}_{\text{BS}}(\phi)$ and $\mathbf{a}_{\text{UE}}(\theta)$ are the normalized array response to an MPC coming from the angles $\phi$ and $\theta$ w.r.t. the BS and UE ULA, respectively. The tuple $(\phi,\theta)\in[-\pi/2,\pi/2]$ is the physical azimuth angle-of-departure (AoD) and angle-of-arrival (AoA), respectively. The normalized ULA responses at the BS and UE are expressed as

[TABLE]

where $d_{\text{BS}}$ and $d_{\text{UE}}$ are the inter-element spacings in the BS and UE ULA, respectively. We assume $d_{\text{BS}}=d_{\text{UE}}=\frac{\lambda}{2}$ where $\lambda$ is the carrier wavelength defined by $\frac{c}{f_{0}}$ with $c$ and $f_{0}$ being the speed of light and the carrier frequency, respectively.

II-C Sparse Beamspace (Virtual) Representation

We use the virtual channel model representation of $\mathbf{H}$, which relates the beamspace and antenna space by the spatial Fourier transform. The virtual channel model describes the channel w.r.t. fixed basis functions corresponding to spatial angles from the finite discrete dictionary. In particular, we follow the framework in Lee et al. [5] and discretize the tuple $(\hat{\phi},\hat{\theta})$ such that the $(\sin(\hat{\phi}),\sin(\hat{\theta}))$ appearing in the array responses (3) are uniformly distributed in [-1,1). Specifically, the quantized grids should satisfy:

[TABLE]

where $G_{\text{UE,D}}$ and $G_{\text{BS,D}}$ are the grid size for the spatial AoA-AoD, respectively. The array response corresponding to the discretized spatial angles are grouped to form the matrices ${\mathbf{A}_{\text{BS,D}}}=\{\mathbf{a}_{\text{BS}}(\hat{\phi});\hat{\phi}\in\mathbf{\Phi}\}$ and ${\mathbf{A}_{\text{UE,D}}}=\{\mathbf{a}_{\text{UE}}(\hat{\theta});\hat{\theta}\in\mathbf{\Theta}\}$, which are transmitting and receiving beamforming matrices at BS and UE, respectively. Generally, $G_{\text{BS,D}}\geq M$ and $G_{\text{UE,D}}\geq N$. However, throughout this work we assume $G_{\text{BS,D}}=M$ and $G_{\text{UE,D}}=N$, resulting in $\mathbf{A}_{\text{BS,D}}\in\mathcal{C}^{M\times M}$ and $\mathbf{A}_{\text{UE,D}}\in\mathcal{C}^{N\times N}$ being unitary Inverse-discrete Fourier transform (DFT) matrices, which are represented as

[TABLE]

The exact beamspace representation can be expressed as $\mathbf{A}_{\text{UE,D}}\mathbf{H}_{\text{V}}\mathbf{A}_{\text{BS,D}}^{\text{H}}$, where $\mathbf{H}_{\text{V}}\in\mathcal{C}^{N\times M}$ is the beamspace sparse matrix defined as follows:

[TABLE]

and $\mathbf{\hat{{a}}}_{\text{UE}}(\theta)$ = $\mathbf{A}_{\text{UE,D}}^{\text{H}}\mathbf{a}_{\text{UE}}(\theta)$ and $\mathbf{\hat{{a}}}_{\text{BS}}(\phi)$ = $\mathbf{A}_{\text{BS,D}}^{\text{H}}\mathbf{a}_{\text{BS}}(\phi)$ are the normalized UE and BS array responses w.r.t. the DFT basis, respectively. The UE array response w.r.t. the DFT basis can be compactly represented as [7, 14],

[TABLE]

where $\vartheta_{n^{\prime},l}=\frac{n^{\prime}-1}{N}-\frac{1}{2}\sin{(\theta_{l})}-\frac{1}{2};n^{\prime}\in[N]$. The proof of the equivalent representation is straightforward and omitted for brevity. Similarly, $\mathbf{\hat{a}}_{\text{BS}}(\phi_{l})$ is formulated as $[\mathbf{\hat{a}}_{\text{BS}}(\phi_{l})]_{m^{\prime}}=\frac{1}{{M}}\frac{\sin(\pi\varphi_{m^{\prime},l}M)}{\sin{(\pi\varphi_{m^{\prime},l})}}e^{-j\pi\varphi_{m^{\prime},l}(M-1)}$, where $\varphi_{m^{\prime},l}=\frac{m^{\prime}-1}{M}-\frac{1}{2}\sin{(\phi_{l})}-\frac{1}{2};m^{\prime}\in[M]$. Substituting $\mathbf{\hat{a}}_{\text{UE}}$ and $\mathbf{\hat{a}}_{\text{BS}}$ into (6), the $(m^{\prime},n^{\prime})$ entry of the beamspace matrix $\mathbf{H}_{\text{V}}$ becomes

[TABLE]

where $\mathcal{D}(\varphi_{m^{\prime},l},\vartheta_{n^{\prime},l})=\frac{1}{{MN}}\frac{\sin(\pi\varphi_{m^{\prime},l}M)}{\sin{(\pi\varphi_{m^{\prime},l})}}\frac{\sin(\pi\vartheta_{m^{\prime},l}N)}{\sin{(\pi\vartheta_{m^{\prime},l})}}$ is the Dirichlet kernel (note that $\mathcal{D}(\varphi_{m^{\prime},l},\vartheta_{n^{\prime},l})=1$ when $\varphi_{m^{\prime},l}=\vartheta_{n^{\prime},l}=0$). Since the DFT is discrete in nature, the beamspace domain $\mathbf{H}_{\text{V}}$ in (8) is evaluated only at integer points ($m^{\prime},n^{\prime}$) with $m^{\prime}\in[M]$, $n^{\prime}\in[N]$. However, the Dirichlet kernel peak need not occur at any of these integer points ($m^{\prime}$, $n^{\prime}$) and can take continuous values, i.e., $m^{\star}\in[1,M],n^{\star}\in[1,N]$. Therefore, the Dirichlet kernel need not peak at one of the pre-defined spatial angles, and maxima of the DFT may not correspond to maxima of the Dirichlet kernel or maxima of the discrete-time Fourier transform (DTFT).

We now sidestep away from the DFT representation to discuss the DTFT concept, which will be pivotal in understanding our algorithms in Section III. The continuum of (8) evaluated at $m^{\star}\in[1,M]$ and $n^{\star}\in[1,N]$ is represented as

[TABLE]

This form (9) is the DTFT counterpart of (8), where $(m^{\star},n^{\star})$ can take continuous values. That is, $1\leq m^{\star}\leq M$ and $1\leq n^{\star}\leq N$. Intuitively, each MPC in the physical domain results in a continuous Dirichlet kernel in the continuum of the beamspace domain, as shown in Fig. 1.

II-D Sparse Recovery Problem

Aided by the sparse virtual representation and vector identity property, vec($\bf{ABC}$) = $(\mathbf{C}^{T}\otimes\mathbf{A})\text{vec}(\mathbf{B})$, MIMO channel estimation (1) is posed as sparse recovery [1] and rewritten:

[TABLE]

where $\mathbf{A}=\mathbf{\Phi\Psi}\in\mathcal{C}^{M_{\text{t}}N_{\text{t}}\times MN}$ is the overall sensing matrix, $\mathbf{\Phi}=(\mathbf{F}^{\text{T}}\otimes\mathbf{W}^{\text{H}})$ represents the combined effect of the precoder and combiner, $\bf{\Psi}=(\mathbf{A}_{\text{BS,D}}^{*}\otimes\mathbf{A}_{\text{UE,D}})$ is the dictionary matrix, and $\bf{A_{\text{BS,D}}}$ and $\mathbf{A}_{\text{UE,D}}$ are matrices for the DFT basis. Finally, $\mathbf{y}\in\mathcal{C}^{M_{\text{t}}N_{\text{t}}\times 1}$ is the vectorized form of $\mathbf{Y}$.

Conventional CS techniques assume that the signal is exactly sparse, which is true only when the physical AoA-AoD tuples are aligned with discretized spatial angles; this is the ideal on-grid case. However, the physical AoA-AoD ($\theta_{l},\phi_{l}$) take continuous values, which may not be aligned with any discretized spatial angles, resulting in basis mismatch (off-grid) effects. These effects violate the sparsity assumption, resulting in performance degradation of CS-based techniques [7]. Below, we analyze the impact of basis mismatch on the sparsity level, which is central to any CS methods.

II-E Effect of Off-Grid/Basis Mismatch

Best Case (on-grid)

The physical AoA-AoD ($\theta_{l},\phi_{l}$) tuple falls exactly on any of the pre-defined spatial angles, and the DFT and DTFT peaks coincide, resulting in exact sparse representation. This phenomenon is illustrated by a 1-dimensional (1D) normalized Dirichlet kernel (Fig. 1, left).

Worst Case (off-grid)

The physical AoA-AoD ($\theta_{l},\phi_{l}$) tuple results in a virtual AoA-AoD ($\varphi_{l},\vartheta_{l}$), which is exactly in between any adjacent pre-defined virtual angles. That is, the resulting virtual AoA-AoD tuple is at distance $(\frac{1}{2M},\frac{1}{2N})$ away from some predefined virtual angle. In this case, the DTFT and DFT peaks do not coincide, which affects not only the two adjacent cells but the entire grid, with amplitude decaying at rate $1/M$ and $1/N$ in the virtual AoA-AOD, respectively, violating the sparsity assumption. Examples of on-grid and worst off-grid effects in a MIMO system appear in Fig. 1 (middle and right panels, respectively).

To provide insights about the worst case of the off-grid problem, we provide a Lemma for MIMO in the presence of a single unit strength MPC, which resembles Lemma 2 in Gao et al. [15]. Our Lemma provides insights about the number of non-zero indices (which we define as the sparsity level) that need to be recovered by the CS methods to capture $\eta$ percent of the power of an MPC in the virtual domain. We concentrate on the worst-case and upper bound the general case, which can be extended to multiple MPCs as well.

Lemma 1

Let $K$ represent the sparsity level in the virtual AoA-AoD domain. Without loss of generality, we assume $K$ to be a multiple of 4. The power captured by the $K$ strongest elements in the DFT domain is given by $\eta=\frac{P_{K}}{P_{T}}$, where

[TABLE]

The ratio $\eta$ is the power captured by $K$ dominant $K$ DFT elements. For the worst and best case scenarios, Fig. 2 shows the number of non-zeros that must be recovered by CS algorithms in order to capture $\eta$ power for a single MPC with $M=N=16$. In other scenarios, the sparsity level is between the blue and red lines. As evident from the figure, in the worst case scenario CS algorithms must recover significantly more non-zeros, which is inefficient. A better way to capture power in the virtual domain requires that we reach the peak of the DTFT spectrum instead of operating on DFT peaks. In other words, we can capture all the power (100%) of each MPC by identifying the peak location of the Dirichlet kernel.

In light of these observations, our main objective is to find the maxima of the Dirichlet kernel peaks (strength/location) instead of recovering all the non-zero elements in the virtual DFT domain. Mathematically, the peak strength and location can be obtained by solving the following optimization,

[TABLE]

[TABLE]

Note that $\mathbf{H}_{\text{V}}$ in our objective function is still evaluated at integer points (${m^{\prime},n^{\prime}}$). This optimization procedure finds the strength and location of the Dirichlet kernel peaks while minimizing the $\ell_{2}$ norm of the residual between the estimated parameters and measurement vector. The constraints mandate that the search is not just over a finite set of angles defined in (4) but over the entire space. Intuitively, finding maxima locations of the Dirichlet kernel (DTFT spectrum) corresponds to estimating the AoA-AoD ($\theta_{l},\phi_{l}$). Similarly, finding the peak strength of the Dirichlet kernel corresponds to estimating the strength of the MPC $\alpha_{l}$.

The above joint optimization problem (12) is non-convex and in general challenging. The non-convexity mainly arises from the Dirichlet structure; Fig.1 exudes it. For brevity, we omit the proof of non-convexity. Below, we propose simpler and efficient schemes to solve (12) using coarse estimates of DFT points obtained from the greedy OMP algorithm.

III Exploiting Dirichlet Kernel Structure

This section begins by investigating the Dirichlet kernel for a single MPC. After addressing a single MPC, we propose an algorithm that accommodates MIMO and multiple MPCs.

To keep our presentation simple, suppose that a single MPC falls off the grid, resulting in a Dirichlet kernel in the virtual AoA domain (Fig. 3). Denote the DFT and Dirichlet kernel peaks by $m^{\prime}$ and $m^{\star}$, respectively. Due to properties of the DFT and DTFT, the Dirichlet kernel peak lies in the range $[m^{\prime}-1,m^{\prime}+1]$. We compute least square (LS) estimates at locations $\{m^{\prime}-1,m^{\prime},m^{\prime}+1\}$ (integer indices adjacent to $m^{\prime}$) under the constraint of Dirichlet structure, resulting in estimates $[\mathbf{H}_{\text{V}}]_{m^{\prime}-1}$, $[\mathbf{H}_{\text{V}}]_{m^{\prime}}$, and $[\mathbf{H}_{\text{V}}]_{m^{\prime}+1}$. The goal is to traverse to the Dirichlet kernel peak and reconstruct it using the LS estimates. This can be achieved in at least two ways.

III-A Dirichlet OMP-Main Lobe (DOMP-MLb)

Without loss of generality, suppose that the Dirichlet peak is in the range $[m^{\prime},m^{\prime}+1]$. In this case, $|[\mathbf{H}_{\text{V}}]_{m^{\prime}+1}|>|[\mathbf{H}_{\text{V}}]_{m^{\prime}-1}|$, implying that the main-lobe (MLb) is in the range $[m^{\prime},m^{\prime}+1]$. Thus, we can ignore the side lobe (SLb) estimate. Based on the MLb estimates, we can find the location $m^{\star}=m^{\prime}+\delta_{m}$, where $\delta_{m}=\frac{1}{2}\min\left(\frac{[\mathbf{H}_{\text{V}}]_{m^{\prime}}}{[\mathbf{H}_{\text{V}}]_{m^{\prime}+1}},\frac{[\mathbf{H}_{\text{V}}]_{m^{\prime}+1}}{[\mathbf{H}_{\text{V}}]_{m^{\prime}}}\right)$ is the deviation from the DFT index, which is obtained by exploiting the uni-modal symmetric, concave property of the main lobe (Fig. 3).

III-B DOMP-Main and Side Lobe (DOMP-MSLb)

Instead of considering just the MLb estimate, the algorithm becomes more robust by not discarding the SLb estimates, but using them to estimate the Dirichlet peak. Based on the main and side lobe (MSLb) estimates, we can find the location $m^{\star}=m^{\prime}+\delta_{m}$, where $\delta_{m}=\scriptsize{\frac{\tan(\frac{\pi}{M})}{\frac{\pi}{M}}\mathbb{R}\left(\frac{[\mathbf{H_{\text{V}}}]_{m^{\prime}-1}-[\mathbf{H}_{\text{V}}]_{m^{\prime}+1}}{2[\mathbf{H}_{\text{V}}]_{m^{\prime}}-[\mathbf{H}_{\text{V}}]_{m^{\prime}-1}-[\mathbf{H}_{\text{V}}]_{m^{\prime}+1}}\right)}$; details in [16].

Extending the previous ideas to a 2D MIMO problem is straightforward. For the 2D problem, suppose that the DFT peaks occur at ($m^{\prime},n^{\prime}$), where $m^{\prime}$ and $n^{\prime}$ refer to indices in the virtual AoA/AoD domains, respectively. Similarly, the Dirichlet kernel peak occurs at ($m^{\star},n^{\star}$). For each MPC in the MIMO case, we need 5 estimates: the DFT peak index ($m^{\prime},n^{\prime}$) and LS estimate ($[\mathbf{H}_{\text{V}}]_{m^{\prime},n^{\prime}}$) and 4 DFT indices and LS estimates around the DFT peak (denoted by $\kappa$ and $[\mathbf{H}_{\text{V}}]_{\kappa}$ in Algorithm 1, respectively). With this information, one can reconstruct the 2D Dirichlet and repeat the procedure recursively for each MPC (denoted by $l$). Further details of our proposed reconstruction method are summarized in Algorithm 1; some remarks are in order.

Remark 1

In some sense, what we describe in Algorithm 1 is a way to identify the true AoA-AoD for each MPC provided the DFT points of the MLb/SLb surrounding each of the Dirichlet kernels. The correct DFT peak points can often be obtained by the computationally tractable OMP projection strategy; steps 2 through 5 in Algorithm 1. Note that the steps of Algorithm 1 after the projection operation (step 2) differ from standard OMP, which iterates between the projection and LS steps without exploiting the structure. Also, the steps for DOMP-MLb/MSLb are the same, except for the update steps enclosed in the if-else statement (steps 6 through 20).

Remark 2

The inherent disadvantage of the DOMP-MLb/MSLb is that Dirichlet kernels in the virtual domain must not overlap, implying that in the physical domain there cannot be closely spaced MPCs. In their current form, DOMP-MLb/MSLb are suitable only for scenarios such as terahertz communication [15] or a single MPC within each spatially separated cone [17]. To overcome this limitation, we next propose a variant of the above method.

III-C DOMP-Local Optimization (DOMP-LO)

The key idea here is to solve (12) in an iterative fashion for each MPC over a localized space around nearby DFT peaks. If the point $(m^{\prime},n^{\prime})$ where the DFT peaks for each MPC are provided, then the search space for the Dirichlet peak can be reduced to the range $([m^{\prime}-1,m^{\prime}+1],[n^{\prime}-1,n^{\prime}+1])$, which turns the joint optimization problem (12) into a convex problem. The localized optimization problem (12) can then be solved for the global optimum for each MPC separately; see step 4 in Algorithm 2. Upon finding the peak location/strength of the Dirichlet kernel, it will be subtracted from the measurements $\mathbf{y}_{\text{res}}$ using step 5. This gets repeated until the stopping criterion is met. Finally, although DOMP-LO overcomes the drawbacks of DOMP-MLb/MSLb by solving the localized optimization problem (its ability to deal with multiple nearby MPCs), it requires more computation.

IV Numerical Results

In this section, the performance of the three proposed channel estimation algorithms are provided, and compared to standard OMP. The performance is evaluated by the normalized mean square error (NMSE) = $\frac{||\mathbf{\hat{H}}-\mathbf{H}||_{F}^{2}}{||\mathbf{H}||_{F}^{2}}$. We consider a scenario with $M=N=G_{\text{UE}}=G_{\text{BS}}=32$ and 3 MPCs. We deliberately place the MPCs away from the grid points. In particular, we place the AoA-AoD at a distance of ($\pm 0.1\zeta/2M,\pm 0.1\zeta/2N$) from the middle of the randomly chosen adjacent grid points, where $\zeta$ is uniformly distributed in the range $[0,1]$. To avoid overlap in the virtual domain, the MPCs are spatially separated by at least $20^{\circ}$ from each other. This scenario helps evaluate the performance of DOMP-MLb/MSLb as noted in Remark 2. However, unlike the DOMP-MLb/MSLb, the DOMP-LO can be used even when the MPCs are closely spaced. All results presented below are averaged over 50 independent trials.

Fig. 4 and Fig. 5 show how exploiting the Dirichlet kernel improves the channel estimation NMSE performance by accounting for off-grid effects, considering different SNRs and different number of measurements. The performance improvement of the DOMP algorithms is mainly due to their ability to combat off-grid effects by traversing Dirichlet kernel peaks.

Among the DOMP methods, the performance of DOMP-MSLb is slightly better than DOMP-MLb, as it is more robust than just considering the estimates from the MLb [16]. The performance gap between standard OMP and DOMP-MLb (and MSLb) is significant, even in the low-SNR and small number of measurements regimes. Moreover, DOMP-LO (Algorithm 2) outperforms other DOMP methods, although this requires more computation for the gradient and the Lagrange multipliers update in solving the localized problem. The detailed computational complexity analysis of the proposed methods is left for future work.

V Conclusion

In this paper, we proposed low-complexity iterative OMP-based algorithms to exploit the implicit Dirichlet structure in the Fourier domain, and thus combat off-grid effects. The proposed DOMP-MLb/MSLb uses the MLb/SLb estimates to traverse to the Dirichlet peaks. On the other hand, the DOMP-LO method solves localized optimization problems to achieve a similar goal. Numerical results show that, in comparison to standard OMP our proposed algorithms achieve lower reconstruction errors in off-grid scenarios for a wide range of SNR and measurement levels.