PARAFAC-Based Multiuser Channel Parameter Estimation for MmWave Massive MIMO Systems over Frequency Selective Fading Channels

Abstract

Channel estimation is crucial in millimeter wave (mmWave) massive multiple-input multiple-output (MIMO) systems, especially with a few training sequences. To solve the problem of uplink channel estimation in mmWave massive MIMO systems, a PARAFAC-based algorithm is proposed for joint estimation of multiuser channels. The orthogonal frequency divisional multiplexing (OFDM) technique is exploited to combat the frequency selective fading channels. In this paper, the received signal at the base station (BS) is formulated as a third-order parallel factor (PARAFAC) tensor, and then a low-complexity algorithm is designed for fast estimation of the factor matrices related to channel parameters, thus leading to joint estimation of multiuser channel parameters via one-dimensional search. Moreover, the Cramér–Rao Bound (CRB) results for multiuser channel parameters are derived for evaluation. Theorical analysis and numerical results reveal that the algorithm performs well with a few training sequences. Compared with existing algorithms, the proposed algorithm has clear advantages both in estimation accuracy and computational complexity.

1. Introduction

As a key technology for 5G wireless communications, millimeter wave (mmWave) massive multiple-input multiple-output (MIMO) has received significant attention in recent years [1,2]. mmWave techniques provide the ability of fast speed data delivery and wide bandwidth [3], making them valued in information ages. The large number of spatial degrees of freedom created by the massive MIMO contribute to high user throughput and high spectral efficiency for extended mobile broadband (eMBB). Furthermore, massive MIMO can be used to efficiently support ultra-reliable low-latency communication (URLLC) and massive machine-type communication (mMTC) in the various types of Internet of Things (IoT) connectivities [4,5,6]. In mmWave massive MIMO systems, large antenna arrays implemented at the base station (BS) and mobile station (MS) are utilized to supply adequate beamforming gains, which offers a compensation for high signal attenuation at mmWave frequencies [7,8,9]. In such a context, the number of radio frequency (RF) chains is much smaller than that of antennas due to lower hardware cost and power consumption. Considering the RF limitations, the hybrid analog/digital precoding scheme [10,11,12,13,14] is widely adopted to provide larger precoding gains. However, the overall channel state information is hard to observe because of massive antennas and hybrid precoding structures. In practice, mmWave channels exhibit the characteristic of broad-band frequency-selective fading (FSF) owing to the large bandwidth and distinct multipath delays [15]. While the whole bandwidth may be subject to FSF, each sub-band undergoes flat fading. When orthogonal frequency division multiplexing (OFDM) is integrated into massive MIMO architectures, mmWave channels are viewed as the aggregation of a number of parallel flat fading channels, which urges us to achieve large multiplexing and diversity gains [16,17,18] with a large number of spatial degrees of freedom. Nevertheless, the overall system performance still hugely depends on the channel estimation, which is a tricky task for the channel subspace sampling limitation [19,20].

In the aspect of channel estimation in massive MIMO systems, several relevant schemes [21,22,23,24,25] have been developed in various scenarios. In general, massive MIMO channels are efficiently parameterized with a series of associated parameters including fading coefficients, angles of arrival or departure (AoAs/AoDs) and delays. Owing to the intrinsic sparsity of mmWave channels, some approaches [10], Refs. [26,27,28,29] based on compressing sensing (CS) techniques convert the issue of channel estimation into the recovery of a line sparse signal. For example, a conventional CS-based algorithm is proposed for channel estimation with much training overhead reduction in a multiuser mmWave system [28]. Further in the work [10], with a novel hierarchical multi-resolution codebook, an adaptive CS-based algorithm is proposed to estimate the channel parameters both in single-path and multipath mmWave environments. According to the CS-based methods, the continuous parameter space is discretized into a finite of grid points [30,31,32], and the estimated angles are assumed to lie on the prespecified grid points. In fact, actual angles may be not in accordance with this assumption, which seriously degrades the channel estimation performance. In addition, some tensor decomposition-based methods stand out in the fields of mmWave channel estimation. The work [33] achieves the joint channel parameter estimation via parallel factor (PARAFAC) analysis in massive MIMO systems. As the improvement of [33], Ref. [34] proposes a more accurate PARAFAC-based estimator in the presence of pilot contamination. For the joint multiuser channel estimation, [35] proposes a layered pilot transmission scheme and decomposes the problem of multiuser channel estimation into a number of problems of single-user channel estimation. In the aforementioned scenarios [33,34,35], simulations demonstrate that the tensor-based algorithms have better estimation performance compared with these CS-based methods. However, most of the existing works concentrate on the narrowband systems while realistic mmWave channels exhibit the broad-band FSF. The work [36] proposes a simultaneous orthogonal matching pursuit (SOMP) scheme and regards the broadband mmWave channel estimation as a multiple measurement vector (MMV) problem with a common support (i.e., the channel support at different frequencies is assumed to be the same). Unlike the proposed scheme in [36], the work [37] estimates the support of the angle-domain channel at some frequencies independently by the orthogonal matching pursuit (OMP) algorithm and combines them into the common support at all frequencies. Moreover, a CANDECOMP/PARAFAC (CP) decomposition-based method for downlink channel estimation is proposed in a mmWave MIMO-OFDM system, where the received signal at individual user is formulated as a low-rank PARAFAC tensor. So far, to the best of our knowledge, there are few tensor decomposition-based methods for joint multiuser mmWave channel estimation, especially over the FSF channels.

In our paper, we consider the broad-band FSF channel and address the problem of multiuser uplink channel estimation in a massive MIMO system, where the hybrid precoding strategy is adopted. Due to the sparsity of mmWave channels, a PARAFAC decomposition-based scheme is proposed for multiuser mmWave channel estimation. In the successive time slots, all the MSs simultaneously transmit respective pilots to a common BS, and the signals received at the BS can be fitted in a low-rank PARAFAC structure.

Table 1. Some reviewed channel estimation works.

Works	System Scenario	Algorithm
[36]	multiuser mmWave massive MIMO-OFDM system	SOMP
[37]	mmWave massive MIMO-OFDM system	OMP
[38]	mmWave massive MIMO-OFDM system	adaptive CS
[10]	mmWave massive MIMO system	CP decomposition-based method
our work	multiuser mmWave massive MIMO-OFDM system	PARAFAC-based channel parameter estimation

presents some reviewed works along their main distinctive characteristics to better position our present work. Unlike the system [38] where individual users are considered for downlink transmission, our proposed algorithm concentrates on the joint estimation of multiuser uplink channels, which leads to wider identifiability and more relaxed requirements for unique PARAFAC decomposition. The main contributions of our work are summarized as follows.

(1) For the sparsity of mmWave FSF channels, the signals transmitted from multiple users can be modeled as a low-rank PARAFAC structure at the BS, which promotes the joint channel estimation in contrast with [38]. It’s worth noting that the massive MIMO channel is characterized with path gains, AoAs/AoDs and delays. This tensor modeling not only takes full advantage of diversities in space, time and frequency dimensions, but also has the benefits of multidimensional structures.

(2) In this context, a PARAFAC decomposition-based algorithm is proposed for accurate estimates of multiuser mmWave channels parameters including path gains, AoAs/AoDs and delays. The proposed algorithm consists of two stages. In the first stage, the accelerated trilinear alternating least squares (ATALS) algorithm is put forward for efficient PARAFAC decomposition. In the second stage, we make further efforts to extract unknown channel parameters of multiple users from factor matrices. The uniqueness of PARAFAC decomposition is also studied, which gives the guidance on the design of precoding and combining matrices. What’s more, our proposed scheme has wider identifiability in terms of unique PARAFAC decomposition.

(3) To provide a benchmark for system evaluation, the Cramér -Rao Bound (CRB) results for channel parameters of multiple users are derived. Simulations reveal that the mean square errors (MSEs) of channel parameters obtained by our algorithm are nearly close to their corresponding CRB results, which indicates our proposed PARAFAC-based algorithm yields accurate channel estimation. While our proposed scheme is put forward at the assumption of uniform linear array (ULA) antennas, it is also suitable for uniform planar array (UPA) antennas.

(4) The computational complexity of the proposed scheme is analyzed. For comparison, we also analyze the computational complexities of some CS-based channel estimation schemes including the SOMP [36], OMP [37] and adaptive CS [10] algorithms. Numerical results demonstrate that our proposed algorithm outperforms its counterparts in estimation accuracy and computational complexity, particularly with a substantial training overhead reduction. Specially, our proposed scheme still performs well, even though the number of users is large.

The rest of this paper is organized as follows. Section 2 first presents the fundamentals of PARAFAC model and formulates the PARAFAC model of multiuser mmWave massive MIMO-OFDM architecture. The uniquess issue of constructed PARARFAC model is also discussed in Section 2. In Section 3, a two-staged PARAFAC-based algorithm is proposed for multiuser mmWave channel estimation. The CRB derivation are given in Section 4. Section 5 presents some numerical results for the system performance analysis. Finally, some conclusions are drawn in Section 6.

Notation: Scalars, vectors, matrices and tensors are denoted by lowercase, boldface lower case, boldface uppercase and calligraphic uppercase letters, e.g., a, $\mathbf{a}$, $\mathbf{A}$ and $\mathcal{A}$, respectively. ${\mathbf{A}}^T$, ${\mathbf{A}}^H$, ${\mathbf{A}}^{*}$ and ${\mathbf{A}}^{†}$ denote the transpose, conjugate transpose, inverse and pseudo-inverse of matrix $\mathbf{A}$, respectively. The operators ⊙, ⊗ and ∘ denote the Khatri-Rao, Kronecker and outer products, respectively. ${\mathbf{I}}_K$ denotes the $K\times K$ identity matrix and ${∥·∥}_F$ denotes the Frobenius norm. The operator $diag\left(\mathbf{a}\right)$ constructs a diagonal matrix with the vector $\mathbf{a}$ and the operator $vec\left(·\right)$ changes a matrix into a column vector by stacking the columns. For convenience, all the abbreviations used in this paper are given in Abbreviations.

2. System Model

As shown in Figure 1, we consider a multiuser mmWave massive MIMO-OFDM system composed of BS and U MSs. The BS is equipped with N antennas and R RF chains, and the u-th MS is equipped with $M_u$ antennas and $R_u=1$ RF chain. Since hybrid precoding structure is adopted in this scenario, the number of RF chains is less than that of antennas, i.e., $R⩽N$, $R_u⩽M_u$. Specially, the number of antennas for each MS is assumed to be equal, i.e., $M_u=M$, $\forall u=1,\cdots ,U$. The number of subcarriers for transmission is assumed to be K.

In the uplink, all the MSs simultaneously transmit respective pilot signals to the common BS via a set of K subcarriers. As mmWave channels exhibit the inherent feature of sparsity, a geometric channel model with finite scatterings is utilized to account for mmWave channels. In terms of the uplink transmission via the k-th subcarrier, the channel of $L_u$ scatterings from the u-th MS to BS can be expressed as ${\mathbf{H}}_{u,k}\in {\mathbb{C}}^{N\times M}$ in the frequency domain

$${\mathbf{H}}_{u,k}={\displaystyle \sum _{l_u=1}^{L_u}}{\mathbf{H}}_{u,l_u}exp\left(-j2\pi f_s{\tau }_{u,l_u}k/K\right),$$

(1)

with the $l_u$-th path

$${\mathbf{H}}_{u,l_u}={\beta }_{u,l_u}{\mathbf{a}}_{BS}\left({\phi }_{u,l_u}\right){\mathbf{a}}_{MS}^T\left({\psi }_{u,l_u}\right)\in {\mathbb{C}}^{N\times M},$$

(2)

for $l_u=1,\cdots ,L_u$, $u=1,\cdots ,U$, where $f_s$ is the sampling rate, $L_u$ denotes the number of propagations and ${\tau }_{u,l_u}$ denotes the time delay of the $l_u$-th path. In this setting of $L_u$ scatterings, a single propagation path ${\mathbf{H}}_{u,l_u}$ is parameterized by a set of channel parameters $\left({\beta }_{u,l_u},{\phi }_{u,l_u},{\psi }_{u,l_u}\right)$, where ${\beta }_{u,l_u}$, ${\phi }_{u,l_u}$ and ${\psi }_{u,l_u}$ denote the complex gain, AoA and AoD of the $l_u$-th path, respectively. The vectors ${\mathbf{a}}_{BS}\left({\phi }_{u,l_u}\right)$ and ${\mathbf{a}}_{MS}\left({\psi }_{u,l_u}\right)$ are the receive and transmit steering vectors, respectively. For the ULA, the steering vectors are given by

$${\mathbf{a}}_{BS}\left({\phi }_{u,l_u}\right)=\frac{1}{\sqrt{N}}{\left[1,e^{j\left(2\pi /\lambda \right)dsin\left({\phi }_{u,l_u}\right)},\cdots ,e^{j\left(N-1\right)\left(2\pi /\lambda \right)dsin\left({\phi }_{u,l_u}\right)}\right]}^T\in {\mathbb{C}}^{N\times 1},$$

(3)

$${\mathbf{a}}_{MS}\left({\psi }_{u,l_u}\right)=\frac{1}{\sqrt{M}}{\left[1,e^{j(2\pi /\lambda )dsin\left({\psi }_{u,l_u}\right)},\cdots ,e^{j(M-1)(2\pi /\lambda )dsin\left({\psi }_{u,l_u}\right)}\right]}^T\in {\mathbb{C}}^{M\times 1},$$

(4)

where $\lambda$ is signal wavelength and d is the antenna spacing. Assume that ${\mathbf{A}}_{BS,u}=\left[{\mathbf{a}}_{BS}\left({\phi }_{u,1}\right),\cdots ,{\mathbf{a}}_{BS}\left({\phi }_{u,L_u}\right)\right]\in {\mathbb{C}}^{N\times L_u}$, ${\mathbf{A}}_{MS,u}=\left[{\mathbf{a}}_{MS}\left({\psi }_{u,1}\right),\cdots ,{\mathbf{a}}_{MS}\left({\psi }_{u,L_u}\right)\right]\in {\mathbb{C}}^{M\times L_u}$, ${\mathit{\beta }}_{u,k}=\left[{\beta }_{u,1,k},\cdots ,{\beta }_{u,L_u,k}\right]\in {\mathbb{C}}^{1\times L_u}$ with the term ${\beta }_{u,l_u,k}={\beta }_{u,l_u}exp\left(-j2\pi f_s{\tau }_{u,l_u}k/K\right)$, we derive

$${\mathbf{H}}_{u,k}={\mathbf{A}}_{BS,u}diag\left({\mathit{\beta }}_{u,k}\right){\mathbf{A}}_{MS,u}^T,$$

(5)

2.1. PARAFAC Model

In this subsection, we give a brief overview of PARAFAC model and discuss the unique issue of PARAFAC decomposition. The PARAFAC decomposition amounts to decomposing the N-way tensor $\mathcal{X}\in {\mathbb{C}}^{I_1\times \cdots \times I_N}$ into a sum of components, each component being a rank-one tensor. For a three-way tensor $\mathcal{X}\in {\mathbb{C}}^{I\times J\times K}$ of R components, PARAFAC decomposition can be expressed as

$$\mathcal{X}={\displaystyle \sum _{r=1}^R}{\mathbf{a}}_r\circ {\mathbf{b}}_r\circ {\mathbf{c}}_r\iff {\left[\mathcal{X}\right]}_{i,j,k}={\displaystyle \sum _{r=1}^R}{\left[\mathbf{A}\right]}_{i,r}{\left[\mathbf{B}\right]}_{j,r}{\left[\mathbf{C}\right]}_{k,r},$$

(6)

where $\mathbf{A}=\left[{\mathbf{a}}_1,\cdots ,{\mathbf{a}}_R\right]$, $\mathbf{B}=\left[{\mathbf{b}}_1,\cdots ,{\mathbf{b}}_R\right]$ and $\mathbf{C}=\left[{\mathbf{c}}_1,\cdots ,{\mathbf{c}}_R\right]$. For brevity, we have the expression $\mathcal{X}=\left[\left[\mathbf{A},\mathbf{B},\mathbf{C}\right]\right]$. Alternatively, it can be written as the mode-n product:

$$\mathcal{X}={\mathcal{I}}_{3,R}{\times }_1\mathbf{A}{\times }_2\mathbf{B}{\times }_3\mathbf{C},$$

(7)

where ${\mathcal{I}}_{3,R}$ is the three-way identity tensor of dimensions $R\times R\times R$ which represents a diagonal hypercubic tensor with its nonzero elements ${\delta }_{i,j,k}=1$ when and only when $i=j=k$. Furthermore, slicing $\mathcal{X}$ along each dimension comes to three distinct slabs $\mathcal{X}\left(i,:,:\right)=\mathbf{B}{D}_i\left(\mathbf{A}\right){\mathbf{C}}^T$, $\mathcal{X}\left(:,j,:\right)=\mathbf{A}{D}_j\left(\mathbf{B}\right){\mathbf{C}}^T$ and $\mathcal{X}\left(:,:,k\right)=\mathbf{A}{D}_k\left(\mathbf{C}\right){\mathbf{B}}^T$. Arranging the slices in some manner leads to the three corresponding matrix unfoldings: ${\mathbf{X}}^{\left(1\right)}=\left(\mathbf{C}\odot \mathbf{B}\right){\mathbf{A}}^T\in {\mathbb{C}}^{KJ\times I}$, ${\mathbf{X}}^{\left(2\right)}=\left(\mathbf{C}\odot \mathbf{A}\right){\mathbf{B}}^T\in {\mathbb{C}}^{IK\times J}$ and ${\mathbf{X}}^{\left(3\right)}=\left(\mathbf{B}\odot \mathbf{A}\right){\mathbf{C}}^T\in {\mathbb{C}}^{IJ\times K}$.

Extended to a N-way tensor $\mathcal{X}\in {\mathbb{C}}^{I_1\times \cdots \times I_N}$ of rank R, PARAFAC analysis has the form:

$$\mathcal{X}={\mathcal{I}}_{N,R}{\times }_{n=1}^N{\mathbf{A}}^{\left(n\right)},$$

(8)

with its element ${\left[\mathcal{X}\right]}_{i_1,\cdots ,i_N}={\displaystyle \sum _{r=1}^R}{\displaystyle \prod _{n=1}^N}{\left[{\mathbf{a}}_r^{\left(n\right)}\right]}_{i_n}$, where ${\mathbf{A}}^{\left(n\right)}\in {\mathbb{C}}^{I_N\times R}$, $i_n=1,\cdots ,I_N$, $n=1,\cdots ,N$, $r=1,\cdots ,R$. Furthermore, the columns of matrix ${\mathbf{A}}^{\left(n\right)}$ are assumed to be normalized, i.e., with a unit norm, for $1⩽n⩽N$, we derive that

$${\left[\mathcal{X}\right]}_{i_1,\cdots ,i_N}={\displaystyle \sum _{r=1}^R}{g}_r{\displaystyle \prod _{n=1}^N}{\left[{\mathbf{A}}^{\left(n\right)}\right]}_{i_n,r},$$

(9)

where the identity tensor ${\mathcal{I}}_{N,R}$ is replaced by the diagonal tensor $\mathcal{G}\in {\mathbb{C}}^{R\times \cdots \times R}$ whose nonzero elements are equal to scaling factors $g_r$ lying on its diagonal.

Then we investigate the uniqueness of PARAFAC-N model in (9). Since the condition for essential uniqueness of PARAFAC model involves the notion of k-rank of a matrix, we first recall the definition of k-rank. The k-rank of a matrix $\mathbf{A}\in {\mathbb{C}}^{I\times R}$ is denoted by $k_{\mathbf{A}}$, where $k_{\mathbf{A}}=r$ on condition that any r columns of $\mathbf{A}$ are linearly independent, where k-rank stands for Kruskal-rank. It can be shown by using the identifiability theorem for the PARAFAC-N model if

$${\displaystyle \sum _{r=1}^R}{k}_{{\mathbf{A}}^{\left(n\right)}}⩾2R+(N-1),$$

(10)

then $\mathcal{X}\in {\mathbb{C}}^{I_1\times \cdots \times I_N}$ is essentially unique and its loading matrices ${\mathbf{A}}^{\left(n\right)}\in {\mathbb{C}}^{I_n\times R}$, $n=1,\cdots ,N$ are unique to column permutation and scaling. Hence, there exists an alternative set of matrices ${\widehat{\mathbf{A}}}^{\left(n\right)}\in {\mathbb{C}}^{I_n\times R}$ with the relation ${\widehat{\mathbf{A}}}^{\left(n\right)}={\mathbf{A}}^{\left(n\right)}\mathbf{\Pi }{\mathbf{\Delta }}^{\left(n\right)}$, for $n=1,\cdots ,N$, where $\mathbf{\Pi }$ is a permutation matrix and ${\mathbf{\Delta }}^{\left(n\right)}$ is a nonsingular diagonal matrix such as ${\displaystyle \prod _{n=1}^N}{\mathbf{\Delta }}^{\left(n\right)}={\mathbf{I}}_R$.

2.2. Constructed PARAFAC Model

According to the upstream communication protocol, multiple MSs encode their respective pilots to the BS through the FSF channels via K subcarriers in P successive time slots. At the destination, the BS employs the combining vectors to detect the transmitted signals from the MSs. Define $s_{u,k}\in {\mathbb{C}}^{R_u\times 1}$ as the pilot for the u-th MS with k-th the subcarrier, and ${\mathbf{e}}_{u,k}^{\left(p\right)}\in {\mathbb{C}}^{M\times R_u}$ as the precoding vector in the p-th time slot, the signal vector received at the BS in the p-th time slot can be expressed as

$${\mathbf{x}}_k^{\left(p\right)}={\displaystyle \sum _{u=1}^U}{\mathbf{H}}_{u,k}{\mathbf{e}}_{u,k}^{\left(p\right)}{s}_{u,k}+{\mathbf{v}}_k^{\left(p\right)},$$

(11)

where ${\mathbf{v}}_k^{\left(p\right)}$ is additive white Gaussian noise. Considering the hybrid precoding strcuture, the BS employs the combining vectors ${\left\{{\mathbf{f}}_{r,k}^{\left(p\right)}\right\}}_{r=1}^R$ to detect the received signal ${\mathbf{x}}_k^{\left(p\right)}$. For the r-th RF chain at the BS, we derive

$$y_{r,k}^{\left(p\right)}={\mathbf{f}}_{r,k}^{{\left(p\right)}^T}{\displaystyle \sum _{u=1}^U}{\mathbf{H}}_{u,k}{\mathbf{e}}_{u,k}^{\left(p\right)}{s}_{u,k}+v_{r,k}^{\left(p\right)},$$

(12)

where $y_{r,k}^{\left(p\right)}$ is the received signal observed from the r-th RF via the k-th subcarrier in the p-th time slot, and $v_{r,k}^{\left(p\right)}$ is additive White Gaussian noise. Arrange the R combining vectors at the BS to form the composite combining matrix ${\mathbf{F}}_k^{\left(p\right)}\in {\mathbb{C}}^{R\times N}$, i.e., ${\mathbf{F}}_k^{\left(p\right)}={\left[{\mathbf{f}}_{1,k}^{\left(p\right)},\cdots ,{\mathbf{f}}_{R,k}^{\left(p\right)}\right]}^T$, which is composed of baseband combiner ${\mathbf{F}}_{BB,k}^{\left(p\right)}\in {\mathbb{C}}^{R\times R}$ and RF combiner ${\mathbf{F}}_{RF}^{\left(p\right)}\in {\mathbb{C}}^{N\times R}$, i.e., ${\mathbf{F}}_k^{\left(p\right)}={\mathbf{F}}_{RF}^{\left(p\right)}{\mathbf{F}}_{BB,k}^{\left(p\right)}\in {\mathbb{C}}^{N\times R}$. Similarly, the combined precoder ${\mathbf{e}}_{u,k}^{\left(p\right)}$ at the u-th MS is composed of baseband precoder $e_{MM,u,k}^{\left(p\right)}\in {\mathbb{C}}^{R_u\times R_u}$ and RF precoder ${\mathbf{e}}_{RF,u}^{\left(p\right)}\in {\mathbb{C}}^{M\times R_u}$, i.e., ${\mathbf{e}}_{u,k}^{\left(p\right)}={\mathbf{e}}_{RF,u}^{\left(p\right)}{e}_{MM,u,k}^{\left(p\right)}\in {\mathbb{C}}^{M\times R_u}$. It should be mentioned that the elements of RF precoding/combining matrices ${\mathbf{e}}_{RF,u}^{\left(p\right)}$ and ${\mathbf{F}}_{RF}^{\left(p\right)}$ meet the constant modulus property. In addition, all K subcarriers share the same RF precoding/combining in consideration of flat analog precoders which are fixed over the band.

Our goal is to achieve reliable estimation of channel parameters $\left({\beta }_{u,l_u},{\phi }_{u,l_u},{\psi }_{u,l_u},{\tau }_{u,l_u}\right)$ from $y_{r,k}^{\left(p\right)}$ with as few pilots as possible. Since the elements of matrices ${\mathbf{e}}_{u,k}^{\left(p\right)}$ and ${\mathbf{F}}_k^{\left(p\right)}$ are randomly chosen from a unit cicrle, we assume that ${\mathbf{e}}_{u,k}^{\left(p\right)}={\mathbf{e}}^{\left(p\right)}$ and ${\mathbf{F}}_k^{\left(p\right)}=\mathbf{F}$. Substituting (5) into (12) yields

$$\begin{array}{cc}\hfill y_{r,k}^{\left(p\right)}& ={\mathbf{f}}_r^T{\displaystyle \sum _{u=1}^U}{\mathbf{A}}_{BS,u}diag\left({\mathit{\beta }}_{u,k}\right){\mathbf{A}}_{MS,u}^T{\mathbf{e}}^{\left(p\right)}{s}_{u,k}+v_{r,k}^{\left(p\right)}\hfill \\ \hfill & ={\mathbf{f}}_r^T{\mathbf{A}}_{BS}diag\left({\mathit{\beta }}_k\right){\mathbf{A}}_{MS}^T{\mathbf{e}}^{\left(p\right)}+v_{r,k}^{\left(p\right)},\hfill \end{array}$$

(13)

where $\Gamma ={\sum }_{u=1}^U{L}_u$, ${\mathbf{A}}_{BS}=\left[{\mathbf{A}}_{BS,1},\cdots ,{\mathbf{A}}_{BS,U}\right]\in {\mathbb{C}}^{N\times \Gamma }$, ${\mathbf{A}}_{MS}=\left[{\mathbf{A}}_{MS,1},\cdots ,{\mathbf{A}}_{MS,U}\right]\in {\mathbb{C}}^{M\times \Gamma }$, ${\mathit{\beta }}_k=\left[{\mathit{\beta }}_{1,k}{s}_{1,k},\cdots ,{\mathit{\beta }}_{U,k}{s}_{U,k}\right]\in {\mathbb{C}}^{1\times \Gamma }$. Define ${\mathbf{Z}}_k={\mathbf{A}}_{BS}diag\left({\mathit{\beta }}_k\right){\mathbf{A}}_{MS}^T$, ${\mathbf{Z}}_k$ corresponds to the frontal slice of third order PARAFAC model $\mathcal{Z}\in {\mathbb{C}}^{N\times M\times K}$. We derive that

$$\mathcal{Z}={\mathcal{I}}_{3,R}{\times }_1{\mathbf{A}}_{BS}{\times }_2{\mathbf{A}}_{MS}{\times }_3\mathbf{B},$$

(14)

where $\mathbf{B}={\left[{{\mathit{\beta }}_1}^T,\cdots ,{{\beta }_K}^T\right]}^T\in {\mathbb{C}}^{K\times \Gamma }$. Clearly, the aggregate channel tensor $\mathcal{Z}\in {\mathbb{C}}^{N\times M\times K}$ includes all the channels of U MSs via K subcarriers with its three modes standing for the number of antennas at the BS, antennas at the MS and subcarriers for transmission.

Collecting the received signal $y_{r,k}^{\left(p\right)}$ from R RF chains contributes to ${\mathbf{y}}_k^{\left(p\right)}={\left[y_{1,k}^{\left(p\right)},\cdots ,y_{R,k}^{\left(p\right)}\right]}^T\in {\mathbb{C}}^{R\times 1}$, and arranging the vectors $\left\{{\mathbf{y}}_k^{\left(1\right)},\cdots ,{\mathbf{y}}_k^{\left(P\right)}\right\}$ in columns leads to the received signal at the k-th subcarrier. By concatenating the received data via K subcarriers, a three-way signal tensor $\mathcal{Y}\in {\mathbb{C}}^{R\times P\times K}$ can be expressed as

$$\mathcal{Y}={\mathcal{I}}_{3,R}{\times }_1{\mathbf{F}}_{BS}{\times }_2{\mathbf{E}}_{MS}{\times }_3\mathbf{B}+\mathcal{V},$$

(15)

where ${\mathbf{F}}_{BS}=\mathbf{F}{\mathbf{A}}_{BS}\in {\mathbb{C}}^{R\times \Gamma }$, $\mathbf{B}\in {\mathbb{C}}^{K\times \Gamma }$, ${\mathbf{E}}_{MS}={\mathbf{E}}^T{\mathbf{A}}_{MS}\in {\mathbb{C}}^{P\times \Gamma }$ and $\mathbf{E}=\left[{\mathbf{e}}^{\left(1\right)},\cdots ,{\mathbf{e}}^{\left(P\right)}\right]\in {\mathbb{C}}^{M\times P}$. Notice that $\mathcal{Y}\in {\mathbb{C}}^{R\times P\times K}$ and $\mathcal{V}\in {\mathbb{C}}^{R\times P\times K}$ are the signal and noise tensors characterized with three dimensions standing for numbers of RF chains, time slots and subcarriers in space, time and frequency domains, respectively. Moreover, the three forms of unfolding ${\mathbf{Y}}^{\left(1\right)}$, ${\mathbf{Y}}^{\left(2\right)}$ and ${\mathbf{Y}}^{\left(3\right)}$ are given by

$$\begin{array}{cc}\hfill {\mathbf{Y}}^{\left(1\right)}& =(\mathbf{B}\odot {\mathbf{E}}_{MS}){\mathbf{F}}_{BS}^T+{\mathbf{V}}^{\left(1\right)}\in {\mathbb{C}}^{PK\times R},\hfill \\ \hfill {\mathbf{Y}}^{\left(2\right)}& =(\mathbf{B}\odot {\mathbf{F}}_{BS}){\mathbf{E}}_{MS}^T+{\mathbf{V}}^{\left(2\right)}\in {\mathbb{C}}^{RK\times P},\hfill \\ \hfill {\mathbf{Y}}^{\left(3\right)}& =({\mathbf{E}}_{MS}\odot {\mathbf{F}}_{BS}){\mathbf{B}}^T+{\mathbf{V}}^{\left(3\right)}\in {\mathbb{C}}^{RP\times K},\hfill \end{array}$$

(16)

where the slices ${\mathbf{V}}^{\left(1\right)}\in {\mathbb{C}}^{PK\times R}$, ${\mathbf{V}}^{\left(2\right)}\in {\mathbb{C}}^{RK\times P}$ and ${\mathbf{V}}^{\left(3\right)}\in {\mathbb{C}}^{RP\times K}$ represent the three forms of unfolding of tensor $\mathcal{V}\in {\mathbb{C}}^{R\times P\times K}$.

Since the number of resolvable paths $K_{\mathrm{path}}$ from MS to BS approximately satisfies the Poisson distribution $K_{\mathrm{path}}\sim max\{Poisson\left(\lambda \right),1\}$ with $\lambda =1.8$ when mmWave works at 28 GHz [39], each MS can be assumed to have 1 or 2 scatterings in the mmWave communications. Thus, the value of $\Gamma$ is usually relatively smaller to the value of each dimension $\left\{R,P,K\right\}$, which enables the signal tensor $\mathcal{Y}\in {\mathbb{C}}^{R\times P\times K}$ to exhibit a low-rank PARAFAC structure.

2.3. Uniqueness Issue

According to the above derivations, it has been proven that the received signal tensor $\mathcal{Y}\in {\mathbb{C}}^{R\times P\times K}$ has a low-rank PARAFAC structure. In this subsection, we study the uniqueness of PARAFAC analysis and provide some sufficient conditions for the unique decomposition of $\mathcal{Y}\in {\mathbb{C}}^{R\times P\times K}$.

Inequality (10) establishes the sufficient condition for the uniqueness of PARAFAC-N model. In the case of our mmWave massive MIMO architecture, if the sum of k-ranks satisfies

$$k_{{\mathbf{F}}_{BS}}+k_{{\mathbf{E}}_{MS}}+k_{\mathbf{B}}⩾2\Gamma +2,$$

(17)

there exists the unique triple $\left({\widehat{\mathbf{F}}}_{BS},{\widehat{\mathbf{E}}}_{MS},\widehat{\mathbf{B}}\right)$ up to permutation and scaling ambiguities, i.e., ${\widehat{\mathbf{F}}}_{BS}={\mathbf{F}}_{BS}\mathbf{\Pi }{\mathbf{\Delta }}^{\left(1\right)}$, ${\widehat{\mathbf{E}}}_{MS}={\mathbf{E}}_{MS}\mathbf{\Pi }{\mathbf{\Delta }}^{\left(2\right)}$, $\widehat{\mathbf{B}}=\mathbf{B}\mathbf{\Pi }{\mathbf{\Delta }}^{\left(3\right)}$, where ${\mathbf{\Delta }}^{\left(1\right)}{\mathbf{\Delta }}^{\left(2\right)}{\mathbf{\Delta }}^{\left(3\right)}=\mathbf{I}$ and ${\widehat{\mathbf{F}}}_{BS}$, ${\widehat{\mathbf{E}}}_{MS}$ and $\widehat{\mathbf{B}}$ are unique estimates of ${\mathbf{F}}_{BS}$, ${\mathbf{E}}_{MS}$ and $\mathbf{B}$. In the generic case, the sufficient condition becomes

$$min\left(R,\Gamma \right)+min\left(P,\Gamma \right)+min\left(K,\Gamma \right)⩾2\Gamma +2.$$

(18)

This indicates that matrices ${\mathbf{F}}_{BS}$, ${\mathbf{E}}_{MS}$ and $\mathbf{B}$ must be of full k-rank for the reason that the k-rank is always less than or equal to the rank, i.e., $k_{\mathbf{A}}⩽rank\left(\mathbf{A}\right)⩽min\left(I,J\right)$, for any matrix $\mathbf{A}\in {\mathbb{C}}^{I\times J}$.

Then we make a profound study of the k-ranks of $\left({\mathbf{F}}_{BS},{\mathbf{E}}_{MS},\mathbf{B}\right)$. Considering the ULA implemented at the BS and MSs, matrices ${\mathbf{A}}_{BS}$ and ${\mathbf{A}}_{MS}$ have the Vandermonde structure, and thus we have $k_{{\mathbf{A}}_{BS}}=min\left(N,\Gamma \right)$, $k_{{\mathbf{A}}_{MS}}=min\left(M,\Gamma \right)$. For the requirement of full k-rank of ${\mathbf{F}}_{BS}$ and ${\mathbf{E}}_{MS}$, an ingenious design of combining matrix $\mathbf{F}$ and precoding matrix $\mathbf{E}$ becomes an indispensable factor. In practice, both ${\mathbf{F}}_{RF}$ and ${\mathbf{e}}_{RF}^{\left(p\right)}$ are accomplished by analog phase shifters, which only change the phase of signals. Therefore, the constant modulus constraint is imposed on the elements of ${\mathbf{F}}_{RF}$ and ${\mathbf{e}}_{RF}^{\left(p\right)}$. Since the total transmission power constraint is enforced by normalizing $e_{MM}^{\left(p\right)}$ such that ${∥{\mathbf{e}}_{RF}^{\left(p\right)}{e}_{MM}^{\left(p\right)}∥}_F^2=M$, we assume that the entries of combining matrix $\mathbf{F}$ and precoding matrix $\mathbf{E}$ are chosen uniformly from a unit circle by a constant, i.e., ${\left[\mathbf{F}\right]}_{i_1,j_1}=\frac{1}{\sqrt{N}}{e}^{j{\upsilon }_{i_1,j_1}}$, $1⩽i_1⩽R$, $1⩽j_1⩽N$, ${\left[\mathbf{E}\right]}_{i_2,j_2}=\frac{1}{\sqrt{M}}{e}^{j{\upsilon }_{i_2,j_2}}$, $1⩽i_2⩽M$, $1⩽j_2⩽P$, where ${\upsilon }_{i_1,j_1}$ and ${\upsilon }_{i_2,j_2}$ follow the independent and identically distribution (i.i.d) uniform distribution $\mathcal{U}\sim \left[0,2\pi \right)$. Thus, the entries of ${\mathbf{F}}_{BS}$ and ${\mathbf{E}}_{MS}$ are viewed as i.i.d Gaussian variables, and we have $k_{{\mathbf{F}}_{BS}}=min\left(R,\Gamma \right)$, $k_{{\mathbf{E}}_{MS}}=min\left(P,\Gamma \right)$. In addition, due to the shift-invariant vector $\mathbf{d}\left({\tau }_{u,l_u}\right)$ of all U MSs, $\mathbf{B}$ is full k-rank, i.e., $k_{\mathbf{B}}=min\left(K,\Gamma \right)$, where $\mathbf{d}\left({\tau }_{u,l_u}\right)={\left[exp\left(-j2\pi f_s{\tau }_{u,l_u}/K\right),\cdots ,exp\left(-j2\pi f_s{\tau }_{u,l_u}K/K\right)\right]}^T$ can be seen as a frequency-domain steering vector pointing towards the delay ${\tau }_{u,l_u}$.

In practice, the number of propagation paths $\Gamma$ is largely limited for mmWave circumstances of poor scatterings. Moreover, the scatterings appear in groups with similar delays, AoAs and AoDs because of massive antennas, which limits effective numbers of active paths. Generally, the number of subcarriers K is larger than $\Gamma$, which turns the uniqueness condition into

$$\begin{array}{c}\hfill \begin{array}{c}min\left(R,\Gamma \right)+min\left(P,\Gamma \right)+min\left(K,\Gamma \right)\hfill \\ =min\left(R,\Gamma \right)+min\left(P,\Gamma \right)+\Gamma ⩾2\Gamma +2\to min\left(R,\Gamma \right)+min\left(P,\Gamma \right)⩾\Gamma +2.\hfill \end{array}\end{array}$$

(19)

From this condition, we can conclude that the BS has the potiential to achieve the acurrate channel estimation in as few as two successive time slots, with the help of $\Gamma$ combining vectors. Hence, our proposed scheme can realize the accurate multiuser channel estimation with much reduced training overhead. However, a slightly larger pair of values for $\left\{R,P\right\}$ is chosen for better performance in our simulations.

3. Proposed PARAFAC Decomposition-Based Channel Estimation Scheme

In the presence of inescapable white noise, we are seeking for some appropriate solutions to the channel estimation issue. To make the most of low-rank PARAFAC strcuture of tesnsor $\mathcal{Y}$, an efficient and low-complexity multiuser channel estimation via PARAFAC decomposition is put forward. The algorithm is divided into two stages. In the first stage, three factor matrices $\left({\widehat{\mathbf{F}}}_{BS},{\widehat{\mathbf{E}}}_{MS},\widehat{\mathbf{B}}\right)$ including channel parameters are uniquely decomposed from $\mathcal{Y}$ by the proposed ATALS algorithm. In the second stage, channel parameters $\left({\widehat{\beta }}_{u,l_u},{\widehat{\phi }}_{u,l_u},{\widehat{\psi }}_{u,l_u},{\widehat{\tau }}_{u,l_u}\right)$ are extracted from these factor matrices through the one-dimensional search method.

3.1. The ATALS Algorithm

As the workhorse technique for low-rank decomposition of three and higher dimensions, the alternating least squares (ALS) algorithm is widely applied in the field of PARAFAC analysis. In order to estimate the channel associated matrices ${\mathbf{F}}_{BS}$, ${\mathbf{E}}_{MS}$ and $\mathbf{B}$ in an efficient way, the ATALS algorithm is proposed as a speed-up version of the classical ALS algorithm. In our simulations, all the channel state information is assumed to be unknown from the MSs or the BS, and the numbers of propagation paths at each MS, the matrices $\mathbf{E}$ and $\mathbf{F}$ are provided in advance.

In terms of the ALS algorithm, the criteria is that when estimating each matrix in the LS algorithm, the other matrices should maintain the previous estimation until the following cost function converges to the minimum value. At the i-th iteration, the cost function ${\delta }^{\left(i\right)}$ is expressed as:

$${\delta }^{\left(i\right)}={∥{\mathbf{Y}}^{\left(1\right)}-({\widehat{\mathbf{B}}}^{\left(i\right)}\odot {\widehat{\mathbf{E}}}_{MS}^{\left(i\right)}){\mathbf{F}}_{BS}^{{\left(i\right)}^T}∥}_F^2,$$

(20)

where ${\widehat{\mathbf{F}}}_{BS}^{\left(i\right)}$, ${\widehat{\mathbf{E}}}_{MS}^{\left(i\right)}$ and ${\widehat{\mathbf{B}}}^{\left(i\right)}$ are the estimated values of ${\mathbf{F}}_{BS}$, ${\mathbf{E}}_{MS}$ and $\mathbf{B}$ acquired from the i-th iteration, respectively. At the i-th iteration, the factor matrices are updated as follows:

$${\mathbf{F}}_{BS}^{\left(i\right)}={\left[{\left({\widehat{\mathbf{B}}}^{(i-1)}\odot {\widehat{\mathbf{E}}}_{MS}^{(i-1)}\right)}^{†}{\mathbf{Y}}^{\left(1\right)}\right]}^T,{\mathbf{E}}_{MS}^{\left(i\right)}={\left[{\left({\widehat{\mathbf{B}}}^{(i-1)}\odot {\widehat{\mathbf{F}}}_{BS}^{\left(i\right)}\right)}^{†}{\mathbf{Y}}^{\left(2\right)}\right]}^T,{\mathbf{B}}^{\left(i\right)}={\left[{\left({\widehat{\mathbf{E}}}_{MS}^{\left(i\right)}\odot {\widehat{\mathbf{F}}}_{BS}^{\left(i\right)}\right)}^{†}{\mathbf{Y}}^{\left(3\right)}\right]}^T$$

(21)

where ${\widehat{\mathbf{F}}}_{BS}^{\left(i-1\right)}$, ${\widehat{\mathbf{E}}}_{MS}^{\left(i-1\right)}$ and ${\widehat{\mathbf{B}}}^{\left(i-1\right)}$ are the estimated values of ${\mathbf{F}}_{BS}$, ${\mathbf{E}}_{MS}$ and $\mathbf{B}$ at the $\left(i-1\right)$-th iteration, respectively. However, in some cases, simulations reveal that the convergence is slow and the estimates of ${\mathbf{F}}_{BS}$, ${\mathbf{E}}_{MS}$ and $\mathbf{B}$ show an increasing tendency towards a certain direction.

To speed up the convergence of the ALS algorithm, a line search process of calculating the linear interpolations for $\left({\mathbf{F}}_{BS},{\mathbf{E}}_{MS},\mathbf{B}\right)$ is carried out before the ALS algorithm. The interpolated values for the i-th iteration are expressed as

$$\begin{array}{cc}\hfill {\widehat{\mathbf{F}}}_{BS}^{\left(new\right)}& ={\widehat{\mathbf{F}}}_{BS}^{\left(i-2\right)}+\rho ({\widehat{\mathbf{F}}}_{BS}^{\left(i-1\right)}-{\widehat{\mathbf{F}}}_{BS}^{\left(i-2\right)}),\hfill \end{array}$$

(22)

$$\begin{array}{cc}\hfill {\widehat{\mathbf{E}}}_{MS}^{\left(new\right)}& ={\widehat{\mathbf{E}}}_{MS}^{\left(i-2\right)}+\rho ({\widehat{\mathbf{E}}}_{MS}^{\left(i-1\right)}-{\widehat{\mathbf{E}}}_{MS}^{\left(i-2\right)}),\hfill \end{array}$$

(23)

$$\begin{array}{cc}\hfill {\widehat{\mathbf{B}}}^{\left(new\right)}& ={\widehat{\mathbf{B}}}^{\left(i-2\right)}+\rho ({\widehat{\mathbf{B}}}^{\left(i-1\right)}-{\widehat{\mathbf{B}}}^{\left(i-2\right)}),\hfill \end{array}$$

(24)

where $\rho$ is the relaxation factor and the differences ${\mathbf{W}}_{{\mathbf{F}}_{BS}}^{\left(i\right)}={\widehat{\mathbf{F}}}_{BS}^{\left(i-1\right)}-{\widehat{\mathbf{F}}}_{BS}^{\left(i-2\right)}$, ${\mathbf{W}}_{{\mathbf{E}}_{MS}}^{\left(i\right)}={\widehat{\mathbf{E}}}_{MS}^{\left(i-1\right)}-{\widehat{\mathbf{E}}}_{MS}^{\left(i-2\right)}$ and ${\mathbf{W}}_{\mathbf{B}}^{\left(i\right)}={\widehat{\mathbf{B}}}^{\left(i-1\right)}-{\widehat{\mathbf{B}}}^{\left(i-2\right)}$ denote the search directions at the i-th iteration. To speed up the convergence, we choose an appropriate step size $\rho (\rho >1)$ [40], which is viewed as $\rho =1$ in the classic ALS algorithm. Assume that $\rho =i^{1/\alpha }$ and $\alpha =3$, the cost function comes to

$$\begin{array}{cc}\hfill {\delta }_{\rho }^{\left(i\right)}& ={∥{\mathbf{Y}}^{\left(1\right)}-({\widehat{\mathbf{B}}}^{\left(new\right)}\odot {\widehat{\mathbf{E}}}_{MS}^{\left(new\right)}){\mathbf{F}}_{BS}^{{\left(new\right)}^T}∥}_F^2\hfill \\ \hfill & ={∥{\mathbf{Y}}^{\left(1\right)}-(({\widehat{\mathbf{B}}}^{\left(i-2\right)}+\rho {\mathbf{W}}_{\mathbf{B}}^{\left(i\right)})\odot ({\widehat{\mathbf{E}}}_{MS}^{\left(i-2\right)}+\rho {\mathbf{W}}_{{\mathbf{E}}_{MS}}^{\left(i\right)})){({\mathbf{F}}_{BS}^{\left(i-2\right)}+\rho {\mathbf{W}}_{{\mathbf{F}}_{BS}}^{\left(i\right)})}^T∥}_F^2.\hfill \end{array}$$

(25)

The relaxation factor $\rho =i^{1/\alpha }$ is accepted on the condition that ${\delta }_{\rho }^{\left(i\right)}⩽{\delta }^{\left(i\right)}$, and then ${\widehat{\mathbf{E}}}_{MS}^{\left(new\right)}$, ${\widehat{\mathbf{E}}}_{MS}^{\left(new\right)}$ and ${\widehat{\mathbf{B}}}^{\left(new\right)}$ are chosen for the i-th iteration in the following ALS algorithm. Otherwise, ${\widehat{\mathbf{F}}}_{BS}^{\left(i-1\right)}$, ${\widehat{\mathbf{E}}}_{MS}^{\left(i-1\right)}$ and ${\widehat{\mathbf{B}}}^{\left(i-1\right)}$ are chosen for further estimates at the i-th iteration. Nevertheless, if this acceleration experiences three consecutive failures, a smaller step size should be chosen as $\rho =i^{1/(\alpha +1)}$.

3.2. Channel Parameter Extraction

On the basis of estimates $\left({\widehat{\mathbf{F}}}_{BS},{\widehat{\mathbf{E}}}_{MS},\widehat{\mathbf{B}}\right)$ from the previous stage, we aim to extract the channel parameters from these estimated factor matrices in this subsection. With the permutation and scaling ambiguities, we have

$${\widehat{\mathbf{F}}}_{BS}={\mathbf{F}}_{BS}\mathbf{\Pi }{\mathbf{\Delta }}^{\left(1\right)},{\widehat{\mathbf{E}}}_{MS}={\mathbf{E}}_{MS}\mathbf{\Pi }{\mathbf{\Delta }}^{\left(2\right)},\widehat{\mathbf{B}}=\mathbf{B}\mathbf{\Pi }{\mathbf{\Delta }}^{\left(3\right)},$$

(26)

where $\mathbf{\Pi }$ is the permutation ambiguity. ${\mathbf{\Delta }}^{\left(1\right)}$, ${\mathbf{\Delta }}^{\left(2\right)}$ and ${\mathbf{\Delta }}^{\left(3\right)}$ are the diagonal scaling matrices with ${\mathbf{\Delta }}^{\left(1\right)}{\mathbf{\Delta }}^{\left(2\right)}{\mathbf{\Delta }}^{\left(3\right)}=\mathbf{I}$. To achieve the estimation of the channel parameters of all U MSs, we resort to a one-dimensional search. For the $l_u$-th path of the u-th MS, we obtain the estimated channel parameters $\left({\widehat{\phi }}_{u,l_u},{\widehat{\psi }}_{u,l_u},{\widehat{\tau }}_{u,l_u}\right)$:

$${\widehat{\phi }}_{u,l_u}=\underset{{\phi }_{u,l_u}}{argmax}\frac{\left|{\left[{\widehat{\mathbf{F}}}_{BS}{\mathbf{e}}_{u,l_u}\right]}^H\left[{\mathbf{a}}_{BS}\left({\phi }_{u,l_u}\right)\right]\right|}{{∥{\widehat{\mathbf{F}}}_{BS}{\mathbf{e}}_{u,l_u}∥}_2{∥{\mathbf{a}}_{BS}\left({\phi }_{u,l_u}\right)∥}_2},$$

(27)

$${\widehat{\psi }}_{u,l_u}=\underset{{\psi }_{u,l_u}}{argmax}\frac{\left|{\left[{\widehat{\mathbf{E}}}_{MS}{\mathbf{e}}_{u,l_u}\right]}^H\left[{\mathbf{a}}_{MS}\left({\psi }_{u,l_u}\right)\right]\right|}{{∥{\widehat{\mathbf{E}}}_{MS}{\mathbf{e}}_{u,l_u}∥}_2{∥{\mathbf{a}}_{MS}\left({\psi }_{u,l_u}\right)∥}_2},$$

(28)

$${\widehat{\tau }}_{u,l_u}=\underset{{\tau }_{u,l_u}}{argmax}\frac{\left|{\left[\widehat{\mathbf{B}}{\mathbf{e}}_{u,l_u}\right]}^H\left[\mathbf{d}\left({\tau }_{u,l_u}\right)\right]\right|}{{∥\widehat{\mathbf{B}}{\mathbf{e}}_{u,l_u}∥}_2{∥\mathbf{d}\left({\tau }_{u,l_u}\right)∥}_2},$$

(29)

where ${\mathbf{e}}_{u,l_u}$ is the l-th unit coordinate vector with $l={\displaystyle \sum _{i=1}^{u-1}}{L}_i+l_u$. Relying on the estimates ${\widehat{\phi }}_{u,l_u}$ and ${\widehat{\psi }}_{u,l_u}$ of U MSs, we can obtain the scaling matrices ${\mathbf{\Delta }}^{\left(1\right)}$ and ${\mathbf{\Delta }}^{\left(2\right)}$. According to the equality ${\mathbf{\Delta }}^{\left(1\right)}{\mathbf{\Delta }}^{\left(2\right)}{\mathbf{\Delta }}^{\left(3\right)}=\mathbf{I}$, matrix ${\mathbf{\Delta }}^{\left(3\right)}$ can be computed. With the known parameters ${\tau }_{u,l_u}$ and ${\mathbf{\Delta }}^{\left(3\right)}$, the complex gain ${\beta }_{u,l_u}$ fianlly can be obtained. So far, we are capable of estimating the channel parameters $\left({\widehat{\beta }}_{u,l_u},{\widehat{\phi }}_{u,l_u},{\widehat{\psi }}_{u,l_u},{\widehat{\tau }}_{u,l_u}\right)$ of U MSs, and then the channel tensor $\mathcal{Z}$ can be recovered. The detailed steps of our proposed algorithm are given in Algorithm 1:

Algorithm 1 PARAFAC decomposition-based channel estimation algorithm

Input: Received tensor $\mathcal{Y}$, precoding matrix $\mathbf{E}$, combining matrix $\mathbf{F}$, one-dimensional search number J, the number of paths $L_u$ of U MSs and error threshold $\epsilon$.

Output: Estimates ${\widehat{\beta }}_{u,l_u}$, ${\widehat{\phi }}_{u,l_u}$, ${\widehat{\psi }}_{u,l_u}$ and ${\widehat{\tau }}_{u,l_u}$.

First stage (the ATALS algorithm):

Initialization: Randomly initialize ${\widehat{\mathbf{F}}}_{BS}^{\left(0\right)}$, ${\widehat{\mathbf{F}}}_{BS}^{\left(1\right)}$, ${\widehat{\mathbf{B}}}^{\left(0\right)}$, ${\widehat{\mathbf{B}}}^{\left(1\right)}$, ${\widehat{\mathbf{E}}}_{MS}^{\left(0\right)}$, ${\widehat{\mathbf{E}}}_{MS}^{\left(1\right)}$. Set $i=2$, $\alpha =3$ and ${\delta }^{\left(0\right)}=\infty$.

Step 1.1. Compute ${\delta }_{\rho }^{\left(i\right)}$ from (25) with $\rho =i^{1/\alpha }$, and compare it with ${\delta }^{\left(i\right)}$ from (20).

If ${\delta }_{\rho }^{\left(i\right)}⩽{\delta }^{\left(i\right)}$, construct ${\widehat{\mathbf{F}}}_{BS}^{\left(new\right)}$, ${\widehat{\mathbf{E}}}_{MS}^{\left(new\right)}$, ${\widehat{\mathbf{B}}}^{\left(new\right)}$ with $\rho =i^{1/\alpha }$; else, construct ${\widehat{\mathbf{F}}}_{BS}^{\left(new\right)}$, ${\widehat{\mathbf{E}}}_{MS}^{\left(new\right)}$, ${\widehat{\mathbf{B}}}^{\left(new\right)}$ with $\rho =1$.

Step 1.2. Update ${\widehat{\mathbf{F}}}_{BS}^{\left(i\right)}$ using ${\widehat{\mathbf{B}}}^{\left(new\right)}$, ${\widehat{\mathbf{E}}}_{MS}^{\left(new\right)}$, ${\mathbf{Y}}^{\left(1\right)}$ by ${\widehat{\mathbf{F}}}_{BS}^{\left(i\right)}={\left[{({\widehat{\mathbf{B}}}^{\left(new\right)}\odot {\widehat{\mathbf{E}}}_{MS}^{\left(new\right)})}^{†}{\mathbf{Y}}^{\left(1\right)}\right]}^T$.

Step 1.3. Update ${\widehat{\mathbf{E}}}_{MS}^{\left(i\right)}$ using ${\widehat{\mathbf{B}}}^{\left(new\right)}$, ${\widehat{\mathbf{F}}}_{BS}^{\left(i\right)}$, ${\mathbf{Y}}^{\left(2\right)}$ by ${\widehat{\mathbf{E}}}_{MS}^{\left(i\right)}={\left[{({\widehat{\mathbf{B}}}^{\left(new\right)}\odot {\widehat{\mathbf{F}}}_{BS}^{\left(i\right)})}^{†}{\mathbf{Y}}^{\left(2\right)}\right]}^T$.

Step 1.4. Update ${\widehat{\mathbf{B}}}^{\left(i\right)}$ using ${\widehat{\mathbf{E}}}_{MS}^{\left(i\right)}$, ${\widehat{\mathbf{F}}}_{BS}^{\left(i\right)}$, ${\mathbf{Y}}^{\left(3\right)}$ by ${\widehat{\mathbf{B}}}^{\left(i\right)}={\left[{({\widehat{\mathbf{E}}}_{MS}^{\left(i\right)}\odot {\widehat{\mathbf{F}}}_{BS}^{\left(i\right)})}^{†}{\mathbf{Y}}^{\left(3\right)}\right]}^T$.

Step 1.5. Compute the new error ${\delta }^{\left(i\right)}$. If $\left|{\delta }^{\left(i\right)}-{\delta }^{\left(i-1\right)}\right|/\phantom{\left|{\delta }^{\left(i\right)}-{\delta }^{\left(i-1\right)}\right|{\delta }^{\left(i\right)}}{\delta }^{\left(i\right)}⩽\epsilon$, then end.

Otherwise, set $i=i+1$ and go to Step 1.1..

Second stage (Channel parameter extraction via one-dimension search):

for$u=1,\cdots ,U$, and for$l_u=1,\cdots ,L_U$

Step 2.1. Estimate the AoA ${\widehat{\phi }}_{u,l_u}$ by (27) with ${\phi }_j\in \left[0,2\pi \right]$, $1⩽j⩽J$.

Step 2.2. Estimate the AoD ${\widehat{\psi }}_{u,l_u}$ by (28) with ${\psi }_j\in \left[0,2\pi \right]$, $1⩽j⩽J$.

Step 2.3. Estimate the delay ${\widehat{\tau }}_{u,l_u}$ by (29) with ${\tau }_j\in \left[0,{\tau }_{max}\right]$, $1⩽j⩽J$.

Step 2.4. Compute the diagonal matrices ${\mathbf{\Delta }}^{\left(1\right)}$, ${\mathbf{\Delta }}^{\left(2\right)}$ and ${\mathbf{\Delta }}^{\left(3\right)}$.

Step 2.5. Estimate the path gain ${\widehat{\beta }}_{u,l_u}$ according to $\widehat{\mathbf{B}}=\mathbf{B}\mathbf{\Pi }{\mathbf{\Delta }}^{\left(3\right)}$.

4. Cramér–Rao Bound Deriavtion

In this section, the derivation of CRB results of channel parameters $\left({\beta }_{u,l_u},{\phi }_{u,l_u},{\psi }_{u,l_u},{\tau }_{u,l_u}\right)$ is given. As known to all, the CRB is a lower bound of the variance of arbitrary unbiased estimator [41], which in turn makes it a reliable benchmark. In this system, the derived CRB results provide a lower bound on the channel prediction error. Therein, the covariance of any biased estimator of the canonical parameter vector $\eta$ satisfies

$$Cov\left(\eta \right)⩾{\mathbf{\Omega }}^{-1}\left(\eta \right),$$

(30)

where $\mathbf{\Omega }\left(\eta \right)$ is the Fisher information matrix.

In view of three equivalent forms of $\mathcal{Y}$ in (16), the entries of noise matrices ${\mathbf{V}}^{\left(1\right)}$, ${\mathbf{V}}^{\left(2\right)}$ and ${\mathbf{V}}^{\left(3\right)}$ are i.i.d zero mean, circularly symmetric complex Gaussian noise, of variance ${\sigma }^2$. As a result, the log-likelihood function of signal tensor $\mathcal{Y}$ is expressed as

$$\begin{array}{cc}\hfill lnL\left(\mathcal{Y}\right)& =-RPKln\left(\pi {\sigma }^2\right)-\frac{1}{{\sigma }^2}{\displaystyle \sum _{r=1}^R}{∥{\mathbf{y}}_r^{\left(1\right)}-(\mathbf{B}\odot {\mathbf{E}}_{MS}){\mathbf{f}}_r^{\left(1\right)}∥}_2^2\hfill \\ \hfill & =-RPKln\left(\pi {\sigma }^2\right)-\frac{1}{{\sigma }^2}{\displaystyle \sum _{p=1}^P}{∥{\mathbf{y}}_p^{\left(2\right)}-(\mathbf{B}\odot {\mathbf{F}}_{BS}){\mathbf{e}}_p^{\left(2\right)}∥}_2^2\hfill \\ \hfill & =-RPKln\left(\pi {\sigma }^2\right)-\frac{1}{{\sigma }^2}{\displaystyle \sum _{k=1}^K}{∥{\mathbf{y}}_k^{\left(3\right)}-({\mathbf{E}}_{MS}\odot {\mathbf{F}}_{BS}){\mathbf{b}}_k^{\left(3\right)}∥}_2^2,\hfill \end{array}$$

(31)

where ${\mathbf{y}}_r^{\left(1\right)}$, ${\mathbf{y}}_p^{\left(2\right)}$ and ${\mathbf{y}}_k^{\left(3\right)}$ represent the r-th, p-th, k-th column of ${\mathbf{Y}}^{\left(1\right)}$, ${\mathbf{Y}}^{\left(2\right)}$ and ${\mathbf{Y}}^{\left(3\right)}$, respectively; ${\mathbf{f}}_r^{\left(1\right)}$, ${\mathbf{e}}_p^{\left(2\right)}$ and ${\mathbf{b}}_k^{\left(3\right)}$ represent the r-th, p-th, k-th column of ${\mathbf{F}}_{BS}^T$, ${\mathbf{E}}_{MS}^T$ and ${\mathbf{B}}^T$, respectively. The matrices ${\mathbf{F}}_{BS}$, ${\mathbf{E}}_{MS}$ and $\mathbf{B}$ contain unknown parameter sets ${\left\{{\mathit{\phi }}_u\right\}}_{u=1}^U$, ${\left\{{\mathit{\psi }}_u\right\}}_{u=1}^U$ and ${\left\{{\mathit{\beta }}_u,{\mathit{\tau }}_u\right\}}_{u=1}^U$, respectively. Define

$${\mathit{\phi }}_u=[{\phi }_{u,1},\cdots ,{\phi }_{u,L_u}],{\mathit{\psi }}_u=[{\psi }_{u,1},\cdots ,{\psi }_{u,L_u}],{\mathit{\beta }}_u=[{\beta }_{u,1},\cdots ,{\beta }_{u,L_u}],{\mathit{\tau }}_u=[{\tau }_{u,1},\cdots ,{\tau }_{u,L_u}],$$

and collect the vectors $\left({\mathit{\beta }}_u,{\mathit{\phi }}_u,{\mathit{\psi }}_u,{\mathit{\tau }}_u\right)$ for all U MSs, we derive the parameter vector $\mathit{\eta }$

$$\eta =\left[{\mathit{\beta }}_1,\cdots ,{\mathit{\beta }}_U,{\mathit{\phi }}_1,\mathbf{\cdots },{\mathit{\phi }}_U,{\mathit{\psi }}_1,\mathbf{\cdots },{\mathit{\psi }}_U,{\mathit{\tau }}_1,\cdots ,{\mathit{\tau }}_U\right]\in {\mathbb{C}}^{1\times 4\Gamma },$$

which consits of all the unknown channel parameters. Furthermore, the log-likelihood function is written as $lnL\left(\mathit{\eta }\right)=lnL\left(\mathcal{Y}\right)$. Thereby, the complex Fisher information matrix is given by

$$\begin{array}{cc}\hfill \mathbf{\Omega }\left(\mathit{\eta }\right)& =E\left\{{\left(\frac{\partial lnL\left(\mathit{\eta }\right)}{\partial \mathit{\eta }}\right)}^H\left(\frac{\partial lnL\left(\mathit{\eta }\right)}{\partial \mathit{\eta }}\right)\right\}\hfill \\ \hfill & =\left[\begin{array}{cccc}{\mathbf{\Omega }}_{\beta \beta }\hfill & {\mathbf{\Omega }}_{\beta \phi }\hfill & {\mathbf{\Omega }}_{\beta \psi }\hfill & {\mathbf{\Omega }}_{\beta \tau }\hfill \\ {\mathbf{\Omega }}_{\phi \beta }\hfill & {\mathbf{\Omega }}_{\phi \phi }\hfill & {\mathbf{\Omega }}_{\phi \psi }\hfill & {\mathbf{\Omega }}_{\phi \tau }\hfill \\ {\mathbf{\Omega }}_{\psi \beta }\hfill & {\mathbf{\Omega }}_{\psi \phi }\hfill & {\mathbf{\Omega }}_{\psi \psi }\hfill & {\mathbf{\Omega }}_{\psi \tau }\hfill \\ {\mathbf{\Omega }}_{\tau \beta }\hfill & {\mathbf{\Omega }}_{\tau \phi }\hfill & {\mathbf{\Omega }}_{\tau \psi }\hfill & {\mathbf{\Omega }}_{\tau \tau }\hfill \end{array}\right]\hfill \end{array}$$

(32)

Any matrix block of $\mathbf{\Omega }\left(\mathit{\eta }\right)\in {\mathbb{C}}^{4\Gamma \times 4\Gamma }$ can be expressed as ${\mathbf{\Omega }}_{xy}\in {\mathbb{C}}^{\Gamma \times \Gamma }$, where x and y represent any symbol of $\left(\beta ,\phi ,\psi ,\tau \right)$. The $(l_1,l_2)$-th element of ${\mathbf{\Omega }}_{xy}$ is expressed as

$${\left[{\mathbf{\Omega }}_{xy}\right]}_{l_1,l_2}=E\left\{\frac{\partial lnL\left(\mathit{\eta }\right)}{\partial x_{u_1,l_{u_1}}^{*}}\frac{\partial lnL\left(\mathit{\eta }\right)}{\partial y_{u_2,l_{u_2}}}\right\}.$$

(33)

where $l_1={\displaystyle \sum _{i=1}^{u_1-1}}{L}_i+l_{u_1}$, $l_{u_1}=1,\cdots ,L_{u_1}$, for $u_1=1,\cdots ,U$; $l_2={\displaystyle \sum _{i=1}^{u_2-1}}{L}_i+l_{u_2}$, $l_{u_2}=1,\cdots ,L_{u_2}$, for $u_2=1,\cdots ,U$.

Further precise derivation can be found in the Appendix A.

5. System Performance Analysis

5.1. Analysis of Computational Complexity

In this subsection, we analyze the computational complexity of our proposed PARAFAC decomposition-based algorithm. Since the traditional CS-based methods are common and effective solutions to the problem of mmWave channel estimation, we compare our method with the SOMP, OMP and adaptive CS-based algorithms.

As shown in

Table 2. Computational complexity of respective algorithms.

Algorithm		Complexity
SOMP		$\mathcal{O}\left(RPK{\Gamma }^2\right)+\mathcal{O}\left(RPMUNK\Gamma \right)$
OMP		$\mathcal{O}\left(RPK{\Gamma }^3\right)+\mathcal{O}\left(RPMUNK\Gamma \right)$
Adaptive CS		$\mathcal{O}\left(SRPK{\Gamma }^2\right)+\mathcal{O}\left(SRPMNK\Gamma \right)$
Proposed	First stage	$\mathcal{O}\left(RPK\Gamma +RPK{\Gamma }^2+{\Gamma }^3\right)$
	Second stage	$\mathcal{O}\left(JK\right)$
	Total	$\mathcal{O}\left(RPK\Gamma +RPK{\Gamma }^2+{\Gamma }^3\right)$

, the overall computational complexity of our proposed algorithm consists of two parts. In the first stage, the major computational task is the ATALS algorithm, which is an iterative algorithm in Step 1.2.–Step 1.4. at each iteration. Take the Step 1.2. for example, to compute the matrix ${\widehat{\mathbf{F}}}_{BS}^{\left(i\right)}$ at the i-th iteration, this step has the complexity in order of $\mathcal{O}\left(RPK\Gamma +RPK{\Gamma }^2+{\Gamma }^3\right)$, which is the same for Step 1.3. and Step 1.4.. In the second stage, the complexity mainly lies in the one-dimensional searching (27) - (29) for channel paremeters of each path. The calculation of delay ${\tau }_{u,l_u}$ has the complexity in order of $\mathcal{O}\left(JK\right)$ for each path, similarly for the calculation of AoAs and AoDs. The total number of paths $\Gamma$ is so small that the main complexity order for the one-dimensional search can be ignored. Finally, the overall computational complexity is represented as $\mathcal{O}\left(RPK\Gamma +RPK{\Gamma }^2+{\Gamma }^3\right)$.

As can be seen from Table 2, the complexity of the SOMP, OMP and adaptive CS-based algorithms can finally be written as $\mathcal{O}\left(RPK{\Gamma }^2\right)+\mathcal{O}\left(RPMUNK\Gamma \right)$, $\mathcal{O}\left(RPK{\Gamma }^3\right)+\mathcal{O}\left(RPMUNK\Gamma \right)$ and $\mathcal{O}\left(SRPK{\Gamma }^2\right)+\mathcal{O}\left(SRPMNK\Gamma \right)$, respectively. S is the number of adaptive stages for each path estimation referred in [10]. Considering the smaller value of $\Gamma$, theoretical analysis shows that our proposed algorithm has a lower computational complexity than the SOMP, OMP and adaptive CS-based algorithms. In our simulations, the average run time of the four algorithms is also recorded for comprehensive comparison.

5.2. Simulation Results

In this subsection, numerical results prove the superiority of our proposed algorithm over other CS-based algorithms through computer emulation. For convenience, system parameter settings in our simulations are represented in

Table 3. Simulation configuration of system parameters.

System Parameter	Configuration
BS antennas N/RF chains R	64/10
MS antennas M/RF chains $R_u$	32/1
Transmission subcarriers K	128
Total number of MSs U/paths $\Gamma$	6/8
Carrier frequency $f_c$	28 GHz
Sampling rate $f_s$	0.25 GHz
AoA/AoD ${\phi }_{u,l_u}$/${\psi }_{u,l_u}$	$\left[0,2\pi \right]$
Path gain ${\beta }_{u,l_u}$	$\mathcal{CN}(0,1)$
Delay ${\tau }_{u,l_u}$	$\left[0,{\tau }_{max}\right]$

. In this context, the BS and the MSs are equipped with the ULA with $N=64$ and $M=32$, respectively. The total number of transmission subcarriers is assumed to be $K=128$. The carrier frequency and the sampling rate are set as $f_c=28GHz$ and $f_s=0.25GHz$, respectively. The separation between antenna elements d is half of signal wavelength $\lambda$. According to the wideband geometric channel model, the mmWave channel parameters are set as follows: ${\beta }_{u,l_u}\sim \mathcal{CN}(0,1)$, ${\phi }_{u,l_u}\sim \left[0,2\pi \right]$, ${\psi }_{u,l_u}\sim \left[0,2\pi \right]$, ${\tau }_{u,l_u}\sim \left[0,{\tau }_{max}\right]$ with the maximum delay ${\tau }_{max}=100ns$. In addition, we assume a reasonably large signal-to-noise ratio (SNR) range of 0dB to 30dB at the MSs with $U=6$, $\Gamma =8$, and the number of scatterings between the BS and each MS is 1 or 2. For representing the estimation accuracy of channel parameters, the MSEs are introduced:

$$\begin{array}{cc}\hfill MSE\left(\mathit{\beta }\right)& ={∥\mathit{\beta }-\widehat{\mathit{\beta }}∥}_2^2,MSE\left(\mathit{\phi }\right)={∥\mathit{\phi }-\widehat{\mathit{\phi }}∥}_2^2,\hfill \\ \hfill MSE\left(\mathit{\psi }\right)& ={∥\mathit{\psi }-\widehat{\mathit{\psi }}∥}_2^2,MSE\left(\mathit{\tau }\right)={∥\mathit{\tau }-\widehat{\mathit{\tau }}∥}_2^2.\hfill \end{array}$$

(34)

where $\mathit{\beta }=[{\mathit{\beta }}_1,\cdots ,{\mathit{\beta }}_U]\in {\mathbb{C}}^{1\times \Gamma }$, $\mathit{\phi }=[{\mathit{\phi }}_1,\cdots ,{\mathit{\phi }}_U]\in {\mathbb{C}}^{1\times \Gamma }$, $\mathit{\psi }=[{\mathit{\psi }}_1,\cdots ,{\mathit{\psi }}_U]\in {\mathbb{C}}^{1\times \Gamma }$, $\mathit{\tau }=[{\mathit{\tau }}_1,\cdots ,{\mathit{\tau }}_U]\in {\mathbb{C}}^{1\times \Gamma }$. The vectors $\widehat{\mathit{\beta }}$, $\widehat{\mathit{\phi }}$, $\widehat{\mathit{\psi }}$ and $\widehat{\mathit{\tau }}$ are the estimates of true channel parameter vectors $\mathit{\beta }$, $\mathit{\phi }$, $\mathit{\psi }$ and $\mathit{\tau }$, respectively. Furthermore, the estimation accuracy of multiuser mmWave channels can also be evaluated by the normalized mean square error (NMSE):

$$NMSE={\displaystyle \sum _{k=1}^K}\frac{{∥{\mathbf{Z}}_k-\widehat{{\mathbf{Z}}_k}∥}_F^2}{{∥{\mathbf{Z}}_k∥}_F^2},$$

(35)

where ${\mathbf{Z}}_k$ represents the multiuser channel matrix via the k-th subcarrier and $\widehat{{\mathbf{Z}}_k}$ is its estimated matrix.

In the first example, we examine the performance of our proposed algorithm at different SNR changes with $P=8$. The Figure 2 depicts the MSEs and CRB results of associated channel parameters versus SNR. From Figure 2, we observe that the MSEs of channel parameters exhibit a significant decline on the increase of SNR. As the benchmark for evaluation, the corresponding CRB results have the similar downward trend and perform approximately as a linear function of SNR. As expected, the MSEs performance of channel parameters are very close to their corresponding CRBs, especially in a higher SNR region. The simulation results indicate that our proposed algorithm enjoys an accurate estimation of multiuser channel parameters.

In the second example, we investigate the MSEs performance of our proposed algorithm as a function of the number of subcarriers K, where we set $SNR=30dB$ and $P=8$. The CRB results of different sets of channel parameters are also represented for comparison. It can be seen from Figure 3 that the proposed algorithm yields accurate estimates even for a small value of K. As the value of K increases, the CRB results of a certain set of channel parameters show a downward trend, and the MSEs performance are close to their corresponding CRBs. We also observe that when the number of subcarriers $K⩽2$, our proposed algorithm fails. In fact, the phenomenon is consistent with our previous analysis on unique PARAFAC decomposition. For the case $P=8⩾\Gamma$ and $R=10⩾\Gamma$, the condition (18) is satisfied only if $K⩾2$.

In the third example, we study the NMSE performance of our proposed algorithm when varying the number of subcarriers K, where we assume $SNR=30dB$. The Figure 4 displays the NMSE performance comparison with the SOMP and OMP algorithms for different time slots $P=8$ and $P=16$. Apparently, with the increasing of K, all the three algorithms have a more and more accurate channel estimation performance. Among these algorithms, the proposed PARAFAC-based estimator presents a clear advantage over the SOMP and OMP algorithms, and it has a much lower error floor. In particular, compare the NMSE performance with different P for a certain algorithm, and we can see that it performs better when $P=16$, especially for the proposed PARAFAC-based and OMP algorithm. Clearly, the increasing number of time slots requires more pilots, which leads to a better estimation performance. However, our goal is to reach a reliable channel estimation with as few measurements as possible. Besides, the NMSE curve obtained by our proposed method accords with previous analysis on unique PARAFAC decomposition. For the case $P=8⩾\Gamma$ and $R=10⩾\Gamma$, when $K⩾2$, the uniqueness condition (18) holds. This indicates that our proposed algorithm performs well even at the minimum number of subcarriers $K=2$.

In the fourth example, we investigate the system performance with different numbers of MSs. The impact of different numbers of MSs on the NMSE performance is shown in Figure 5, where the numbers of MSs are assumed to be $U=2$, $U=4$, $U=8$, $U=16$ and $U=32$, with single path propagation for per MS. To guarantee the uniqueness of PARAFAC analysis, we set $P=22$ and $R=12$. As the number of MSs increases, the NMSE performance decreases. There is a reason that channel parameters to be estimated grow with the rise of U, thus leading to a decline in estimation accuracy. Moreover, the estimation performance of factor matrices in the first stage suffers from degradation for the increasing unknown matrix elements. However, our proposed scheme still has a good performance even though the number of MSs is large to $U=32$. While our proposed algorithm is put forward for multiuser architecture, it is also available in a single-user environment of multipaths.

In the final example, we compare our proposed algorithm with the SOMP, OMP and adaptive CS algorithms with different numbers of RF chains $R=8$ and $R=16$. As shown in Figure 6, it describes the NMSE comparison among the four algorithms versus SNR. As the SNR increases, all the algorithms have more and more accurate channel estimation. Our proposed algorithm is superior to the other three CS-based algorithms in terms of estimation accuracy with the same R. Especially in a higher SNR region, the three CS-based algorithms have a poor performance with a higher error floor while the PARAFAC-based method still yields more accurate channel estimates. Such a performance gap comes from that these CS-based algorithms assume the estimated angles are discretely distributed on the pre-defined grids. In fact, actual angles are continuous and may not necessarily lie on the grids, which results in an inevitable grid mismatch, and thus degrading the performance. Generally speaking, the finer the grid is, the smaller the grid mismatch error is. However, it also brings an unstable performance with higher computational complexity. By contrast, our proposed method is not only free of grid mismatch but also fully utilizes the low-rank structure of mmWave channels, which leads to a better and stable performance. From Figure 6, we can also gain some insights into how different numbers of RF chains effect the system performance. As R increases from 8 to 16, the NMSE performance degrades. While more RF chains brings about a better system performance, it suffers from higher hardware cost and power consumption. Therefore, a smaller number of RF chains should be selected for a compromise between hardware cost and system performance.

For this example, the average run time of the SOMP, OMP, adaptive CS and proposed PARAFAC-based algorithms at different SNR changes is shown in

Table 4. Average run time of respective algorithms.

SNR(dB)	0	5	10	15	20	25	30
SOMP (R = 8)	1.1239 s	1.1892 s	1.1923 s	1.2035 s	1.2315 s	1.2174 s	1.2248 s
OMP (R = 8)	2.2667 s	2.2075 s	2.2738 s	2.2493 s	2.2512 s	2.2555 s	2.2542 s
Adaptive CS (R = 8)	1.0299 s	0.9888 s	1.1948 s	1.1877 s	1.2464 s	1.2162 s	1.1911 s
Proposed (R = 8)	0.7117 s	0.6612 s	0.595 s	0.7766 s	0.7377 s	0.6304 s	0.6987 s
SOMP (R = 16)	1.635 s	1.5626 s	1.6055 s	1.6113 s	1.5785 s	1.6153 s	1.6129 s
OMP (R = 16)	5.3028 s	5.3577 s	5.2848 s	5.2799 s	5.328 s	5.4576 s	5.3256 s
Adaptive CS (R = 16)	1.4869 s	1.679 s	1.6707 s	1.678 s	1.6739 s	1.7011 s	1.6774 s
Proposed (R = 16)	0.7749 s	0.8524 s	0.8869 s	0.8478 s	0.8686 s	0.8329 s	0.8868 s

. Numerical results demonstrate that the average run time of our proposed algorithm is shorter than that of other three algorithms. In contrast with the CS-based algorithms, our PARAFAC-based algorithm takes advantage of its inherent structure of multidimensional data and achieves a stable and low-complexity channel estimation.

6. Conclusions

In this paper, a PARAFAC decomposition-based algorithm has been provided for multiuser uplink channel estimation over the broad-band FSF mmWave channels. Owing to the sparsity of mmWave channels, the received signal at the BS is formulated as a low-rank PARAFAC tensor, and thus multiuser channel parameter estimation can be achieved via PARAFAC decomposition. Therein, the uniqueness issue is studied. Moreover, the CRB results are derived for evaluation in the situation of multiple users. Simulations reveal that the algorithm performs well with a few training sequences. Compared with other CS-based algorithms, our proposed algorithm outperforms its counterparts both in estimation accuracy and computational complexity. For future research, we will pay more attention to the time-varying mmWave channels, and try to seek more precise tensor models, such as the nested PARAFAC model and the nested Tucker model. Novel tensor fitting algorithms such as the Bayesian inference-based algorithm [42] will be developed in future work for better performance.