Beam-Delay Domain Channel Estimation for mmWave XL-MIMO Systems

At millimeter-wave frequencies, the channels are typically contributed by very few paths due to the sparse scattering nature [12, 13, 14], which facilitates compressed sensing techniques to enhance the channel estimation performance. For conventional mMIMO systems, the orthogonal matching pursuit (OMP)-based channel estimation algorithm is proposed in [15], which exploits beam domain sparsity based on the angle sampling and reconstructs the channel by extracting the dominant paths. For orthogonal frequency division multiplexing (OFDM) systems, [16] and [17] present the simultaneously weighted OMP (SW-OMP) algorithm and subcarrier selection-simultaneous iterative gridless weighted-orthogonal least squares (SS-SIGW-OLS) to exploit the joint sparsity of beam domain channels under different subcarriers. Furthermore, channel estimation based on beam-delay domain sparsity is proposed in [18] to explore the inter-antenna and inter-subcarrier correlations.

Since the near-field effect in XL-MIMO systems destroys the plane wave assumption, the beam domain sparsity based on the angle sampling no longer holds, necessitating the novel beam domain representation. Following the spherical wave in the near-field channels, the beam domain sparsity based on the angle-distance sampling, denoted as polar domain, is proposed in [19, 20], followed by the OMP algorithm for channel estimation. To alleviate the estimation error introduced by discrete sampling, [19] presents the polar domain simultaneously iterative gridless weighted (P-SIGW) algorithm, which updates the angle-distance sampling grids based on the gradient descent. As the scatterers and MTs may be distributed in both far-field and near-field regions, the hybrid-field model is employed in [21, 22, 23, 24] for such scenarios, which exploit both angle and angle-distance sampling.

With the transmission bandwidth and array aperture growth, the beam-squint effect becomes more significant, which destroys the joint sparsity of beam domain across subcarriers. To address this issue, [25] and [26] propose the pattern of beam domain sparsity variations with subcarriers in far-field channels, followed by the generalized simultaneous OMP (GSOMP) the beam-split pattern detection (BSPD) algorithms. Furthermore, by exploiting the beam-delay domain sparsity with beam-squint effects, the super resolution-based channel estimation algorithms are proposed in [27, 28, 29]. Most recently, the channel estimation with both near-field and beam-squint effects is investigated in [30, 31]. Specifically, [30] explores the bilinear patterns of beam domain sparsity variations with subcarriers in the common sensing matrix, facilitating the bilinear patterns detection (BPD)-based channel estimation algorithm. Besides, the near-field beam-split-aware OMP (NBA-OMP) algorithms is proposed in the [31] based on the GSOMP algorithm, which exploits the frequency selective sensing matrices to explore beam domain sparsity based on angle-distance domain sampling and the joint sparsity across subcarriers. The exploration of joint sparsity variations due to beam squint effects in [31, 30] provides one solution for the channel estimation in mmWave XL-MIMO systems, paving the way for subsequent work.

In mmWave XL-MIMO systems, the near-field and beam-squint effects induce substantial changes in channel models, necessitating more accurate modeling. Besides, there are significant correlations between subcarriers in the OFDM systems, which are under-explored in such scenarios. Finally, the sparse model based on discretized sampling encounters an inherent mismatch with continuous parameters, thereby warranting further investigation into effective mitigation strategies for this issue.

Motivated by the previous works, we investigate the channel estimation for mmWave XL-MIMO systems with near-field and beam-squint effects. The main contributions of this paper are summarized as follows:

• •

  By incorporating the near-field and beam-squint effects, we model the channel in the beam-delay domain to explore inter-antenna and inter-subcarrier correlations. In doing so, we end up with the beam-delay domain sparse representation of the spatial-frequency domain channel. To capture the beam-delay domain sparsity, the independent and non-identically distributed Bernoulli-Gaussian models with unknown prior hyperparameters are employed, facilitating the formulation of the channel estimation problem.

• •

  To cope with the unknown hyperparameters in the prior models, we introduce the constrained Bethe free energy (BFE) minimization framework, which enables the joint estimation of the beam-delay domain channel and the hyperparameters in the prior model. The hybrid message passing (HMP) algorithm is established by designing different belief structures to estimate beam-delay domain channels and their prior hyperparameters, thereby, spatial-frequency domain channels.

• •

  To further improve the model accuracy, we introduce the multidimensional grid point perturbation (MDGPP)-based representation, which assigns individual perturbation parameters to each multidimensional discrete grid and thus mitigates mismatches between continuous parameters and discrete grids. Under the constrained BFE minimization framework, the MDGPP parameters are treated as unknown hyperparameters, facilitating the HMP algorithm for MDGPP-based channel estimation. To reduce the computational complexity, we propose the two-stage HMP algorithm, where the output of the initial estimation stage is pruned for the refinement stage.

Organization: The rest of this paper is organized as below. Section II describes the system model and formulates the channel estimation problem. In Section III, the HMP algorithm is developed based on constrained BFE minimization framework. The two-stage HMP algorithm for the MDGPP-based channel estimation is proposed in Section IV, whose computational complexity is also analyzed. Numerical simulations in Section V shows the superiority of the proposed algorithms over the benchmarks. Finally, Section VI concludes the paper.

Notation: Throughout this paper, imaginary unit is represented by $j=\sqrt{-1}$. Lowercase, bold lowercase, and bold uppercase letters denote scalars, column vectors, and matrices, respectively. The transpose, conjugate, and conjugate-transpose operations are represented by the superscripts $(\cdot)^{T}$, $(\cdot)^{*}$, and $(\cdot)^{H}$, respectively. The Kronecker and Hardmard products are represented by the operator $\otimes$ and $\circ$, respectively. $[\cdot]_{i}$ and $[\cdot]_{i,j}$ are the $i$-th element of a vector and the $(i,j)$-th element of a matrix, respectively. $\norm{\cdot}_{2}$ and $\norm{\cdot}_{\infty}$ represent the $\ell_{2}$ and $\ell_{\infty}$ norms, respectively. $\max\left\{\cdot\right\}$ and $\min\left\{\cdot\right\}$ denote the maximum and minimum operators, respectively. $\mathsf{E}\left\{\cdot\right\}$ and $\mathsf{V}\left\{\cdot\right\}$ denote the expectation and variance operators, respectively. $\lfloor\cdot\rfloor$, $\lceil\cdot\rceil$, and $\langle\cdot\rangle$ denote the floor, ceil, and modulo operations, respectively. The element-wise inverse is represented by the superscripts $(\cdot)^{{\circ}-1}$. The symbols $\mathbb{C}$ denote the complex number fields. $\mathcal{CN}(x;{\mu},{\sigma}^{2})$ denotes the random variable $x$ obeying circularly symmetric complex Gaussian distribution with mean $\mu$ and variance ${\sigma}^{2}$.

In the mmWave XL-MIMO systems, the BS usually adopts hybrid precoding architecture, which makes only a far smaller number of baseband received signals at each training symbol observed. With the hybrid precoding architecture, we assume that the BS is equipped with $N_{\text{RF}}$ RF chains ($N_{\text{RF}}{\;\ll\;}N$), where the RF chains and the antennas are connected by phase shifters network[32]. Besides, quasi-static channel assumptions are considered in this paper, where BS, MTs and scatters remain nearly unchanged during channel estimation. Thus, the received signal $\mathbf{y}_{b,p}$ under the $p$-th pilot subcarrier at the $b$-th training symbol can be expressed as¹¹1Since non-overlapping pilot subcarriers are allocated for different MTs, we can focus on only the single MT and omit the MT index to simplify the expression.

where $x_{b,p}$ denotes the pilot under the $p$-th pilot subcarrier at the $b$-th training symbol, which can be set to $1$ without loss of generality, and $\mathbf{z}_{b,p}{\;\sim\;}\mathcal{CN}\left(0,{\sigma}_{z}^{2}\mathbf{I}\right)$ denotes the additive white Gaussian nosie (AWGN) under the $p$-th pilot subcarrier.

In the hybrid precoding architecture, the number of measurements required for reliable channel estimation cannot be satisfied with only one pilot symbol, and thus $N_{\text{PS}}$ pilot symbols are required. Stacking all the received signal under the $p$-th pilot subcarrier in a column vector, we can obtain

where $\mathbf{F}_{p}=\left[\mathbf{F}_{0,p},{\dots},\mathbf{F}_{N_{\text{RX}}-1,p}\right]{\;\in\;}\mathbb{C}^{N{\times}M}$ denotes the hybrid precoder and $M=N_{\text{RF}}N_{\text{PS}}$ denotes the number of measurements. Then, all received signals under $N_{\text{p}}$ pilot subcarriers are stacked in a column vector, so that the received signal $\mathbf{y}$ can be expressed as

where $\bar{\mathbf{F}}=\mathsf{blkdiag}\{\mathbf{F}_{0},{\dots},\mathbf{F}_{N_{\text{p}}-1}\}$ denotes the hybrid precoder of all pilot subcarriers, and $\mathbf{h}$ denotes the spatial-frequency domain channel.

We assume that the channel between the target MT and the BS has $L$ paths. Then the channel between the MT and the $n$-th antenna of BS under the $p$-th pilot subcarrier can be expressed as

where ${\alpha}_{l}$ denotes the complex gain of $l$-th path, $r_{l}^{\left(n\right)}$ denotes the distance between the $n$-th antenna of BS and the last-hop scatter (LHS)²²2For LOS path and NLOS path, the LHS is defined as the MT and the last scatterer before arriving at the antenna array, respectively. of $l$-th path, $\bar{\tau}_{l}$ denotes the propagation delay between the MT and LHS of $l$-th path, $f_{p}=f_{\text{c}}+\left(p-N_{\text{p}}/2\right){\Delta}\bar{f}$ denotes the operating frequency of the $p$-th pilot subcarrier, and $c$ denotes the speed of light.

As shown in Fig. 1, the distance between the $n$-th antenna of BS and LHS of the $l$-th path for MT satisfies [7]:

where ${\theta}_{l}^{\left(n\right)}$ denotes the physical angle between the $n$-th antenna of ULA and LHS of the $l$-th path, $n_{o}$ denotes the index of reference antenna, ${\Delta}n{\;\triangleq\;}n-n_{\text{o}}$ denotes the antenna index difference between the $n$-th antenna and the reference antenna, and $\left(a\right)$ follows from the fact that the geometric constraint $r_{l}^{\left(n\right)}\cos{\theta}_{l}^{\left(n\right)}$ is constant, $\forall n,l$. With the Taylor series expansion of $\sqrt{1+x}$ at $x=0$, $r_{l}^{\left(n\right)}$ can be approximated as

where ${\psi}_{l}^{\left(n_{\text{o}}\right)}{\;\triangleq\;}\sin\psi_{l}^{\left(n_{\text{o}}\right)}$ denotes the angle between the $n$-th antenna of BS and LHS of the $l$-th path, $d=c/(2f_{\text{c}})$ denotes the half-wavelength antenna spacing, and $\left(a\right)$ is obtained by ignoring the higher order terms of $\Delta{n}d$.

For convenience, $r_{l}{\;\triangleq\;}r_{l}^{\left(n_{\text{o}}\right)}$ and ${\psi}_{l}{\;\triangleq\;}{\psi}_{l}^{\left(n_{\text{o}}\right)}$ are regarded as the location parameters of the $l$-th path for MT without causing confusion. Substituting (II-A) into (4) and stacking the spatial domain channels into vector form, we can obtain

where ${\tau}_{l}=\bar{\tau}_{l}+r_{l}/c$ denotes the propagation delay between the MT and reference antenna at BS, $\mathbf{a}\left(r,{\psi},f_{p}\right)$ denotes phase differences under the $p$-th pilot subcarrier, which can be defined by

$\mathbf{b}\left(r,{\psi}\right)$ denotes power differences, which can be defined by

where $\varrho$ denotes the normalized factor, such that

and the common phase and gain between all antennas independent of ${\Delta}n$ is merged into the path gain $\bar{\alpha}_{l}$. For convenience, we define $\mathbf{c}\left(r,{\psi},f_{p}\right){\;\triangleq\;}\mathbf{a}\left({r},{\psi},f_{p}\right){\;\circ\;}\mathbf{b}\left({r},{\psi}\right)$ as beam domain steering vector.

In (7), the differences between the mmWave XL-MIMO and mMIMO channel model are two-fold:

1. 1.

   Near-field effect: In mmWave XL-MIMO systems, the MTs and LHSs located within the Rayleigh distance result in the spherical wave propagation [6, 7]. Moreover, the propagation distance difference causes fluctuations in power attenuation across antennas.

2. 2.

   Beam-squint effect: In mmWave XL-MIMO systems, the propagation delay differences on the antenna array are not negligible due to the large bandwidth of mmWave systems, which leads to the frequency-dependent beam domain steering vector [8].

To illustrate the deviation from the plane wave assumption due to near-field effects and facilitate subsequent beam domain representations, we can introduce the slope parameter $\eta{\;\triangleq\;}\left(1-{\psi}^{2}\right)/r$ to characterize the rate of instantaneous angle change with ${\Delta}n$, which is shown by

Then, the angle-distance sampling can be converted to the angle-slope sampling by replace the $\left(1-{\psi}^{2}\right)/r$ with ${\eta}$, resulting in the beam domain steering vector $\mathbf{c}\left({\eta},{\psi},f_{p}\right)$.

By stacking the spatial channel of all pilot subcarriers, the spatial-frequency channel is expressed as

where $\mathbf{d}\left({\tau}\right)$ denotes the delay domain steering vector defined as $\left[\mathbf{d}\left({\tau}\right)\right]_{p}=\exp(-j2{\pi}f_{p}{\tau})$, and $\bar{\mathbf{c}}\left({\eta},{\psi}\right)$ denotes aggregate beam domain steering vector given by

From (12), the spatial-frequency domain channel is superposed by several paths, resulting in the inter-antenna and inter-subcarrier correlations. To exploit such correlations, we develop the beam-delay domain sparse representation. Toward this end, we assume that the angles, slopes and delays for each MT satisfy ${\psi}{\;\in\;}\left[-1,1\right)$, ${\eta}{\;\in\;}\left[0,{\eta}_{\max}\right]$, and ${\tau}{\;\in\;}\left[0,{\tau}_{\max}\right]$, respectively. The maximum slope parameter ${\eta}_{\max}=1/r_{\min}$, where $r_{\min}$ denotes the minimum service distance, and ${\tau}_{\max}$ denotes the maximum allowable delay, which is determined by the cyclic prefix length and the array aperture under the beam-squint effect. Then, the angle, slope, and delay domain are sampled uniformly at intervals ${\psi}_{\Delta}=2/N$, ${\eta}_{\Delta}$³³3The sampling interval ${\eta}_{\Delta}$ is determined by the discrete Fresnel function with the coherence threshold [19]., and ${\tau}_{\Delta}=1/(N_{\text{p}}{\Delta}\bar{f})$ to obtain the sampling sets, which can be defined by (14a), (14b), and (14c), respectively.

where $\mathcal{K}_{\text{x}}{\;\triangleq\;}\left\{0,1,{\dots},K_{\text{x}}-1\right\},\text{x}{\;\in\;}\left\{\text{an},\text{sl},\text{de}\right\}$ denotes the sampling point index set, $K_{\text{an}}=\lceil 2/{\psi}_{\Delta}\rceil$, $K_{\text{an}}=\lceil{\eta}_{\max}/{\eta}_{\Delta}\rceil$, and $K_{\text{de}}=\lceil{\tau}_{\max}/{\tau}_{\Delta}\rceil$ denote the number of sampling points in angle, slope, and delay domains, respectively. Since the sampling intervals decrease as $N$ and $N_{\text{p}}$ increase, the spatial-frequency domain channel is well approximated as a linear combination of basis functions determined by angle-slope-delay tuples $\left(\bar{\psi}_{k_{\text{an}}},\bar{\eta}_{k_{\text{sl}}},\bar{\tau}_{k_{\text{de}}}\right)$ for sufficiently large $N$ and $N_{\text{p}}$.

Therefore, with the pre-defined sampling set of the angles, slopes, and delays, (12) can be approximated as

where $\mathbf{u}\left(\bar{\psi}_{k_{\text{an}}},\bar{\eta}_{k_{\text{sl}}},\bar{\tau}_{k_{\text{de}}}\right)$ denotes the angle-slope-delay domain basis function defined as

and ${\beta}_{k_{\text{an}},k_{\text{sl}},k_{\text{de}}}$ denotes the complex gain of the path with respect to angle-slope-delay tuple $\left(\bar{\psi}_{k_{\text{an}}},\bar{\eta}_{k_{\text{sl}}},\bar{\tau}_{k_{\text{de}}}\right)$, which is given by

Reformulating (15) into matrix form and assuming that the spatial-frequency domain channels can be represented perfectly, we can obtain

where $\mathbf{U}{\;\in\;}\mathbb{C}^{NN_{\text{p}}{\times}K_{\text{an}}K_{\text{sl}}K_{\text{de}}}$ and $\bm{\beta}{\;\in\;}\mathbb{C}^{K_{\text{an}}K_{\text{sl}}K_{\text{de}}{\times}1}$ denote the beam-delay domain transformation matrix and the beam-delay domain channel defined as

and

respectively. Without causing confusion, $(k_{\text{an}},k_{\text{sl}},k_{\text{de}})$ maps to the simple index $k{\;\in\;}\{0,1,{\dots},K-1\}$ in a one-to-one manner, which satisfies

and $K=K_{\text{an}}K_{\text{sl}}K_{\text{de}}$.

With the proposed sparse representation, the channel estimation problem is to estimate the beam-delay domain channel and ultimately the spatial-frequency domain channel based on the received signal. The MMSE channel estimator can be expressed as the a posteriori mean [33], given by

where $\mathbf{s}{\;\triangleq\;}\mathbf{G}\bm{\beta}$ and $\mathbf{G}{\;\triangleq\;}\bar{\mathbf{F}}^{H}\mathbf{U}$ denote the auxiliary vector and measurement matrix, respectively, $\mathsf{Pr}\left(\bm{\beta},\mathbf{s}\mid\mathbf{y}\right)$ denotes the posterior probability density function (PDF) of $\bm{\beta}$ and $\mathbf{s}$.

From the received signal model and the beam-delay domain channel representation, $p\left(\bm{\beta},\mathbf{s}\mid\mathbf{y}\right)$ can be expressed as

where $\mathsf{Pr}\left(\mathbf{y}\mid\mathbf{s}\right)$, $\mathsf{Pr}\left(\mathbf{s}\mid\bm{\beta}\right)$, and $\mathsf{Pr}\left(\bm{\beta}\right)$ denote received noise model, baseband transfer model, and beam-delay domain channel prior model, respectively.

Since AWGN channels are employed, the conditional PDF $\mathsf{Pr}\left(\mathbf{y}\mid\mathbf{s}\right)$ can be detailed as

where $y_{m,p}$ and $s_{m,p}$ denote the $\left(mN_{\text{p}}+p\right)$-th element of $\mathbf{y}$ and $\mathbf{s}$, respectively. Then the conditional PDF $\mathsf{Pr}\left(\mathbf{s}\mid\bm{\beta}\right)$ can be factorized as

where $\mathbf{g}_{m,p}^{(\text{r})}$ denotes the $\left(mN_{\text{p}}+p\right)$-th row of $\mathbf{G}$.

With the assumption of wide-sense stationary uncorrelated scattering Rayleigh fading channel, the elements of $\bm{\beta}$ are statistically independent for sufficiently large $N$ and $N_{\text{p}}$ [34]. Moreover, the sparse nature of mmWave channels is considered, implying that only a few number of elements in $\bm{\beta}$ have significant power, with most elements being zero. Thus, we can model $\bm{\beta}$ as an independent and non-identically distributed Bernoulli-Gaussian distribution, and the prior PDF $\mathsf{Pr}\left(\bm{\beta};\bm{\omega}\right)$ controlled by hyperparameters $\bm{\omega}$ can be decomposed into

where $\bm{\omega}=\left[{\lambda},{\chi}_{0},{\dots},{\chi}_{K-1}\right]$, ${\lambda}$ denotes the sparsity level, and ${\chi}_{k}$ denotes the power of ${\beta}_{k}$ when ${\beta}_{k}$ is nonzero.

Building on the probabilistic modeling, the MMSE estimation of spatial-frequency domain channel can be given by

due to the commutability of the MMSE estimator over linear transformations [33].

Following the VFE minimization framework [35, 36, 37], we approximate $\mathsf{Pr}\left(\bm{\beta},\mathbf{s}\mid\mathbf{y};\bm{\omega}\right)$ by the trial belief $\mathsf{b}\left(\bm{\beta},\mathbf{s};\bm{\omega}\right)$, which can be obtained by minimizing the Kullback-Leibler (KL) divergence of $b\left(\bm{\beta},\mathbf{s};\bm{\omega}\right)$ and $\mathsf{Pr}\left(\bm{\beta},\mathbf{s}\mid\mathbf{y};\bm{\omega}\right)$. The KL divergence minimization problem can be expressed by

where $\mathcal{Q}$ denotes the constrained set of $\mathsf{b}\left(\bm{\beta},\mathbf{s};\bm{\omega}\right)$, $D\left[\cdot\|\cdot\right]$ denotes the KL divergence, $\mathcal{F}_{V}$ denotes the VFE defined by

where $\mathsf{Pr}\left(\bm{\beta},\mathbf{y};\bm{\omega}\right)$ denotes the joint PDF, and $-\ln Z\left(\mathbf{y}\right){\;\triangleq\;}-\ln p\left(\mathbf{y}\right)$ denotes the Helmholtz free energy independent of $\bm{\beta}$. Thus, the KL divergence minimization is equivalent to the VFE minimization.

The VFE minimization is still intractable in large-scale systems, which motivates the constrained set design of $\mathsf{b}\left(\bm{\beta};\bm{\omega}\right)$. To balance the exactness and tractability, the Bethe method is adopted to design the constrained set and convert the VFE minimization into the BFE minimization. Following Bayes rule and the probabilistic model, the joint PDF can be factorized as

where $f_{y,m,p}$, $f_{s,m,p}$ and $f_{\beta,k}$ denote the factor nodes. Based on the Bethe method [38, 39, 40], we introduce auxiliary beliefs $\mathsf{b}_{s,m,p}\left(s_{m,p}\right)$, $\mathsf{b}_{s\bm{\beta},m,p}\left(s_{m,p},\bm{\beta}\right)$ and $\mathsf{b}_{{\beta}{\lambda}{\chi},k}\left({\beta}_{k},{\lambda},{\chi}_{k}\right)$ for factor nodes $f_{y,m,p}$, $f_{s,m,p}$ and $f_{{\beta},k}$, respectively. Similarly, we introduce auxiliary beliefs $\mathsf{q}_{\lambda}\left({\lambda}\right)$, $\mathsf{q}_{{\chi},k}\left({\chi}_{k}\right)$, $\mathsf{q}_{{\beta},k}\left({\beta}_{k}\right)$ and $\mathsf{q}_{s,m,p}\left(s_{m,p}\right)$ for variable nodes ${\lambda}$, ${\chi}_{k}$, ${\beta}_{k}$ and $s_{m,p}$, respectively. Following this, the trial belief and the BFE can be expressed as

and

respectively, where $H\left[\cdot\right]$ denotes the entropy.

Inspired by the variational message-passing algorithm [36], we can factorize complicated beliefs into the product of simple beliefs by exploiting the slow-varying characteristic of hyperparameters, which can be expressed as

where $\mathsf{b}_{\beta,k}\left({\beta}_{k}\right)$, $\mathsf{b}_{\lambda}\left({\lambda}\right)$ and $\mathsf{b}_{\chi,k}\left({\chi}_{k}\right)$ denote the auxiliary belief with respect to ${\beta}_{k}$, ${\lambda}$ and ${\chi}_{k}$, respectively. Furthermore, the hyperparameter beliefs can be constrained by the Dirac delta function, which is given as

where ${\lambda}_{\text{true}}{\;\in\;}\left(0,1\right)$ and ${\chi}_{k,\text{true}}{\;\geq\;}0$ denote the ground-truth hyperparameters.

To guarantee the global dependence of factor nodes, the auxiliary beliefs of factor and variable nodes in Bethe method should satisfy marginal consistency constraints (MCCs), which is intractable for continuous random variables. To address this issue, the MCCs can be relaxed to mean and variance consistency constraints (MVCCs), which can be given by

where the beliefs without index subscripts denote the vector version of the scalar beliefs. The variance consistency constraints with respect to $\bm{\beta}$ and $\mathsf{b}_{s\bm{\beta}}$ are further relaxed as (36) by averaging, which can be given by

Finally, the KL divergence minimization problem can be converted into the constrained BFE minimization problem, which is expressed by

The optimization problem in (37) can be solved by the Lagrange multiplier method [41], which sets the derivative of the Lagrangian function with respect to the auxiliary beliefs to zero. Since the auxiliary beliefs of the hyperparameters can be completely factorized in (33), their MCCs can be satisfied constantly and thus can be ignored in the Lagrangian function. Substituting constraints (33) and (34) into $\mathsf{b}_{{\beta}{\lambda}{\chi},k}$, the Lagrangian function of (37) can be expressed as

where $\mathcal{F}_{\text{B, HP}}$ denotes the BFE with factorized hyperparameter constraints, which can be obtained by replacing the KL-divergence $D[\mathsf{b}_{{\beta}{\lambda}{\chi},k}\|f_{\beta,k}]$ with $D[\mathsf{b}_{{\beta},k}\|\hat{f}_{\beta,k}]$ in (III-A), and $\hat{f}_{\beta,k}{\;\triangleq\;}\mathsf{Pr}({\beta}_{k};{\lambda}_{\text{true}},{\chi}_{k,\text{true}})$, $\mathcal{L}_{\text{M}}$ and $\mathcal{L}_{\text{V}}$ defined in (39) denote the subfunctions of mean consistency constraints and variance consistency constraints, respectively, $\bm{\xi}^{s,b_{s}}$, $\bm{\xi}^{s,b_{s\bm{\beta}}}$, $\bm{\xi}^{{\beta},b_{{\beta}}}$, and $\bm{\xi}_{m,p}^{{\beta},b_{s\bm{\beta}}}$ denote the Lagrangian multiplier for mean consistency constraints, $\bm{\varsigma}^{s,b_{s}}$, $\bm{\varsigma}^{s,b_{s\bm{\beta}}}$, $\bm{\varsigma}^{{\beta},b_{{\beta}}}$ and $\bm{\varsigma}^{{\beta},b_{{\beta}}}$ denote the Lagrangian multiplier for variance consistency constraints.

Due to space limitations, the stationary-point equations are provided in Appendix A. Building on the derived equations, the resulting HMP algorithm for channel estimation is summarized in Algorithm 1, where the adaptive damping based on the BFE evaluation [42] is introduced to improve the convergence.

To eliminate the gap between the sampling sets and the continuous parameters, the spatial-frequency domain channel can be decorrelated by the MDGPP-based channel basis functions, which is expressed as

where ${\Delta}{\psi}_{k}{\;\in\;}\left[-{\psi}_{\Delta}/2,{\psi}_{\Delta}/2\right)$, ${\Delta}{\eta}_{k}{\;\in\;}\left[-{\eta}_{\Delta}/2,{\eta}_{\Delta}/2\right)$, and ${\Delta}{\tau}_{k}{\;\in\;}\left[-{\tau}_{\Delta}/2,{\tau}_{\Delta}/2\right)$ denote the perturbation parameters in angle, slope, and delay of the $k$-th sampling tuple, respectively, ${\beta}_{k}$ denotes the complex gain of the $k$-th sampling tuple, and $\tilde{\mathbf{u}}_{k}\left({\Delta}{\psi}_{k},{\Delta}{\eta}_{k},{\Delta}{\tau}_{k}\right)$ denote the perturbed basis function of the $k$-th sampling tuple defined as

where the index tuple $\left({k}_{\text{an}},{k}_{\text{sl}},{k}_{\text{de}}\right)$ and the simple index $k$ satisfies the mapping rule (21). The proposed MDGPP-based representation assigns individual perturbation parameters to each beam-delay domain sampling grid, while the conventional perturbation assigns common perturbation parameters based on grid line constraints [47]. To visualize advantages of the MDGPP-based representation, the schematic of the two-dimensional example is shown in Fig. 2, where the perturbations are shown on the $(\theta,\phi)$ parameter. It can be observed that, with conventional perturbations, the perturbation of “Grid 1” for matching “Path 1” causes the mismatch between “Grid 2” and “Path 2” due to the common perturbation parameter ${\Delta}{\theta}_{n}$, which can be avoided by assigning different perturbation parameters ${\Delta}{\theta}_{n,p}$ and ${\Delta}{\theta}_{n,q}$ in the MDGPP.