{'fraction_non_alphanumeric': 0.046527347490564974, 'fraction_numerical': 0.009772996113332958, 'mean_word_length': 4.198682284040996, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 1, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 1, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
MIMO-OFDM-Based Massive Connectivity With Frequency Selectivity Compensation

In this paper, we study how to efficiently and reliably detect active devices and estimate their channels in a multiple-input multiple-output (MIMO) orthogonal frequency-division multiplexing (OFDM) based grant-free non-orthogonal multiple access (NOMA) system to enable massive machine-type communications (mMTC). First, by exploiting the correlation of the channel frequency responses in narrow-band mMTC, we propose a block-wise linear channel model. Specifically, the continuous OFDM subcarriers in the narrow-band are divided into several sub-blocks and a linear function with only two variables (mean and slope) is used to approximate the frequency-selective channel in each sub-block. This significantly reduces the number of variables to be determined in channel estimation and the sub-block number can be adjusted to reliably compensate the channel frequency-selectivity. Second, we formulate the joint active device detection and channel estimation in the block-wise linear system as a Bayesian inference problem. By exploiting the block-sparsity of the channel matrix, we develop an efficient turbo message passing (Turbo-MP) algorithm to resolve the Bayesian inference problem with near-linear complexity. We further incorporate machine learning approaches into Turbo-MP to learn unknown prior parameters. Numerical results demonstrate the superior performance of the proposed algorithm over state-of-the-art algorithms.


I. INTRODUCTION
M ASSIVE machine-type communications (mMTC) has been envisioned as one of the three key application scenarios of fifth-generation (5G) wireless communications. To support massive connectivity of machine-type devices, mMTC typically has the feature of sporadic transmission with short packets, which is different from conventional humantype communications [1]. This implies that only a small subgroup of devices are active in any time instance of mMTC. As such, in addition to channel estimation and signal detection, a fundamentally new challenge for the design of an mMTC receiver is to reliably and efficiently identify which subgroup of devices are actively engaged in packet transmission.
Recently, a new random access protocol termed grant-free non-orthogonal multiple access (NOMA) has been evaluated and highlighted for mMTC [2]. In specific, in grant-free NOMA, the devices share the same time and frequency resources for signal transmission, and the signals can be transmitted without the scheduling grant from the base station (BS). The receiver at the BS is then required to perform active device detection (ADD), channel estimation (CE), and signal detection (SD), either separately or jointly. The earlier work [3]- [6] assumed that full channel state information (CSI) can be acquired at the BS and studied joint CE and SD. However, the assumption of full CSI availability is not practical since it will cause a huge overhead to estimate the CSI of all access devices. The follow-up work [7] proposed to divide the process at the BS into two phases, namely, the joint ADD and CE phase the SD phase. Since the BS only needs to estimate the CSI of active devices, the pilot overhead is significantly reduced. In addition, the channel sparsity in the device domain enables the employment of compressed sensing (CS) algorithms [8] to solve the joint ADD and CE problem. For example, the authors in [9] considered the multiple-input multiple-output (MIMO) transmission and leveraged a multiple measurement vector (MMV) CS technique termed vector approximate message passing (Vector AMP) [10] to achieve asymptotically perfect ADD. In [11], the authors further considered a mixed analogto-digital converters (ADCs) architecture at the BS antennas and proposed a CS algorithm based on the turbo compressive sensing (Turbo-CS) [12]. It is known from [12] that Turbo-CS outperforms approximate message passing (AMP) [13] both in convergence performance and computational complexity. Another line of research considered the more challenging joint ADD, CE, and SD problem [14], [15]. As compared to the two-phase approach, these joint schemes can achieve significant performance improvement but at the expense of higher computational complexity.
Orthogonal frequency division multiplexing (OFDM) is a mature and enabling technology for 5G to provide high spectral efficiency. As such, the design of OFDM-based grantfree NOMA has attracted much research interest in recent years [16]- [18]. In [16], the authors exploited the blocksparsity of the channel responses on OFDM subcarriers to design a message-passing-based iterative algorithm. Besides, it has been demonstrated in [17] that the message-passingbased iterative algorithm can be unfolded into a deep neural network. By training the parameters of the neural network, the convergence and performance of the algorithm are improved. Furthermore, OFDM-based grant-free NOMA with massive MIMO has been considered in [18]. By leveraging the sparsity both in the device domain and the virtual angular domain, the authors utilized the AMP algorithm to achieve the joint ADD and CE. One issue in [16]- [18] is that the frequencydomain channel estimation on every subcarrier requires a high pilot overhead, which is inefficient for the short data packets transmission in mMTC.
To reduce the pilot overhead, a common strategy is to transform the frequency-domain channel into the time-domain channel by inverse discrete Fourier transform (IDFT) [19]. Due to the limited delay spread in the time-domain channel, the arXiv:2104.05169v1 [cs.IT] 12 Apr 2021 time-domain channel is sparse, thereby requiring fewer pilots for CE. Furthermore, by exploiting the sparsity of both the time-domain channel and the device activity pattern, some state-of-the-art CS algorithms such as Turbo-CS [12], [20] and Vector AMP [9], [10] can be applied to the considered systems with some straightforward modifications. However, there exists an energy leakage problem caused by the IDFT transform to obtain the time-domain channel. The energy leakage compromises the channel sparsity in the time domain. In addition, the power delay profile (PDP) is generally difficult for the BS to acquire, and thus cannot be exploited as prior information to improve the system performance.
Motivated by the bottleneck of the existing channel models when applied to the MIMO-OFDM-based grant-free NOMA system, we aim to construct a channel model to enable efficient massive random access. Due to the short packets in mMTC, the bandwidth for packet transmission is usually narrow, e.g., 1MHz for 10 5 access devices [21]. Then the variations of the channel frequency responses across the subcarriers are limited and slow. By leveraging this fact, we propose a block-wise linear channel model. Specifically, the continuous subcarriers are divided into several sub-blocks. In each subblock, the frequency-selective channel is approximated by a linear function. We demonstrate that the number of variables in the block-wise linear channel model is typically much less than the number of non-zero delay taps in the time-domain channel. Moreover, the sub-block number can be modified to strike the trade-off between the model accuracy and the number of the channel variables to be estimated.
Based on the block-wise linear system model, we aim to design a CS algorithm to solve the joint ADD and CE problem. We first introduce a probability model to characterize the block-sparsity of the channel matrix. Then the joint ADD and CE is formulated as a Bayesian inference problem. Inspired by the success of Turbo-CS [12], [20] in sparse signal recovery, we design a message passing algorithm termed turbo message passing (Turbo-MP) to resolve the Bayesian inference problem. The message passing processes are derived based on the principle of Turbo-CS and the sum-product rule [22]. By designing a partial orthogonal pilot matrix with fast transform, Turbo-MP achieves near-linear complexity.
Furthermore, we show that machine learning methods can be incorporated into Turbo-MP to learn the unknown prior parameters. We adopt the expectation maximization (EM) algorithm [23] to learn the prior parameters. We then show how to unfold Turbo-MP into a neural network (NN), where the prior parameters are seen as the learnable parameters of the neural network. Numerical results show that NN-based Turbo-MP has a faster convergence rate than EM-based Turbo-MP. More importantly, we show that Turbo-MP designed for the propose frequency-domain block-wise linear model significantly outperforms the state-of-the-art counterparts, especially those message passing algorithms designed for the time-domain sparse channel model [9], [11].

A. Notation and Organization
We use bold capital letters like X for matrices and bold lowercase letters x for vectors. (·) T and (·) H are used to denote the transpose and the conjugate transpose, respectively. We use diag(x) for the diagonal matrix created from vector x, diag(x 1 , ..., x N ) for the block diagonal matrix with the n-th block being vector x n , and diag(X 1 , ..., X N ) for the block diagonal matrix with the n-th block being matrix X n . We use vec(X) for the vectorization of matrix X and ⊗ for the Kronecker product. ||X|| F and ||x|| 2 are used to denote the Frobenius norm of matrix X and the l 2 norm of vector x, respectively. Matrix I denotes the identity matrix with an appropriate size. For a random vector x, we denote its probability density function (pdf) by p(x). δ(·) denotes the Dirac delta function and δ[·] denotes the Kronecker delta function. The pdf of a complex Gaussian random vector x ∈ C N with mean m and covariance C is denoted by The remainder of this paper is organized as follows. In Section II, we introduce the existing system models. Furthermore, we propose the block-wise linear system model and demonstrate its superiority. In Section III, we formulate a Bayesian inference problem to address the ADD and CE problem. In Section IV, we propose the Turbo-MP algorithm, describe the pilot design, and analyze the algorithm complexity. In Section V, we apply the EM algorithm to learn the prior parameters. Moreover, we show how Turbo-MP is unfolded into a neural network. In Section VI, we present the numerical results. In Section VII, we conclude this paper.

A. MIMO-OFDM-Based Grant-free NOMA Model
Consider a MIMO-OFDM-based grant-free NOMA system as shown in Fig. 1, where a frequency band of N adjacent OFDM subcarriers are allocated for mMTC. N is usually much less than the total number of available subcarriers in the considered OFDM system. This frequency allocation strategy follows the idea of narrow-band internet-of-things [24] and network slicing [25]. Then the allocated bandwidth is used to support K single-antenna devices to randomly access an Mantenna base station (BS). In each time instance, only a small subset of devices are active. To characterize such sporadic transmission, the device activity is described by an indicator function α k as with p(α k = 1) = λ where λ 1. We adopt a multipath block-fading channel, i.e., the multipath channel response remain constant within the coherence time. Denote the channel frequency response on the n-th subcarrier from the k-th device at m-th BS antenna by where ∆f is the OFDM subcarrier spacing; L k is the number of channel taps of the k-th device; ρ k,l and τ k,l 1 are respectively the l-th tap power and tap delay of the k-th device; β k,m,l ∼ CN (β k,m,l ; 0, 1) is the normalized complex gain and assumed to be independent for any k, m, l [26]. Then the channel frequency response can be expressed in a matrix form as k,n be the pilot symbol of the k-th device transmitted on th n-th subcarrier at the t-th OFDM symbol, and T be the number of OFDM symbols for pilot transmission. Then we construct a block diagonal matrix Λ k ∈ C T N ×N with the n-th diagonal block being [a Assume the cyclic-prefix (CP) length L cp > τ k,l , ∀k, l. After removing the CP and applying the discrete Fourier transform (DFT), the system model in the frequency domain is described as where Y ∈ C T N ×M is the observation matrix; N is an additive white Gaussian noise (AWGN) matrix with its elements independently drawn from CN (0, σ 2 N ). In [16]- [18], CS algorithms were proposed based on (5) to achieve the CE on every OFDM subcarrier.
Define the time-domain channel matrix of the k-th device as H k ∈ C N ×M . Note thatH k can be represented as the IDFT of G k , i.e,H where F ∈ C N ×N is the DFT matrix with the (n 1 , n 2 )-th element being 1/  represented as whereH = [H T 1 , ...,H T K ] T is the channel matrix in the time domain. It is known that the channel delay spread is usually much smaller than N , i.e., some rows ofH k are zeros. Besides, due to the sporadic transmission of the devices, H k is an all-zero matrix if device k is inactive. In this case, CS algorithms such as Vector AMP [9], [10] and Turbo-CS [12], [20] can be used to recoverH from observation Y by exploiting the sparsity ofH. However, path delay is generally not an integer multiple of the sampling interval, resulting in the energy leakage problem of the IDFT transform (6) which severely compromises the sparsity ofH. An illustration of the energy leakage problem is given in Fig. 2(b).

B. Proposed Block-Wise Linear System Model
For narrow-band mMTC [1], the variations of the channel frequency response across the subcarriers are typically slow and limited, which inspires us to develop an alternative channel model to efficiently leverage the correlation in the frequency domain. In specific, we propose a block-wise linear channel model. We divide the N continuous subcarriers into Q subblocks. In each sub-block q, a linear function is used as the approximation of the channel frequency response: where h k,q,m and c k,q,m represent the mean and slope of the linear function in the q-th sub-block, respectively; ∆ k,nq,m is the error term due to model mismatch; and l q = (q− 1 2 )N/Q is the midpoint of n q . Intuitively, h k,q,m can be seen as the meanvalue of the channel response in the q-th sub-block, and c k,q,m is used to characterize the change of the channel response for the compensation of the frequency-selectivity. For the k-th device, we define the matrix H k ∈ C Q×M and C k ∈ C Q×M as The reason for introducing α k in (9)-(10) is that the channel estimation and device activity detection can be jointly achieved by recovering H k and C k .
Then the block-wise linear approximation of the channel frequency response G k is given by where ∆ k ∈ C N ×M is the error matrix from the k-th device with the (n, m)-th element ∆ k,n,m in (8). Substituting (11) into (5) with some manipulations, we obtain the block-wise linear system model as .., C T K ] T ∈ C QK×M is the channel compensation matrix; W is the summation of the AWGN and the error terms from model mismatch as Following the central limit theorem, for a large devices number K, W can be modeled as an AWGN matrix with mean zero and variance σ 2 w . We note that with Q = N and C = 0, system model (12) reduces to the model (5) in [16]- [18] where the channel response on every subcarrier needs to be estimated exactly. Furthermore, we adjust the sub-block number Q to strike a balance between the number of channel variables to be estimated and the model accuracy. An example is shown in Fig. 2 where the channel response in the frequency domain and its block-wise linear approximation is given in Fig. 2(a). It is seen that sub-block number Q = 4 is sufficient to ensure that the block-wise linear model approximates the frequencydomain channel response very well. This implies that, with model (12), only 2Q variables need to be estimated for each device at each BS antenna. Compared to the channel response in the time domain as shown in Fig. 2(b), it is clear that the number of unknown variables 2Q is much less than the number of corresponding non-zero taps. In the remainder of the paper, we focus on the estimation algorithm design based on model (12). In the simulation section, we will show that our algorithm designed based on (12) significantly outperforms the existing approaches based on (5) and (7).

III. PROBLEM STATEMENT
With model (12), our goal is to recover the channel mean matrix H and channel compensation matrix C based on the noisy observation Y. This task can be constructed as a Bayesian inference problem. In the following, we first introduce the probability model of H and C, and then describe the statistical inference problem.
Due to the sporadic transmission of the devices, matrices H and C have a structured sparsity referred to as block-sparsity. In specific, if the k-th device is inactive, we have H k = 0 and C k = 0. With some abuse of notation, we utilize a conditional Bernoulli-Gaussian (BG) distribution [19] to characterize the block-sparsity as where δ(·) is the Dirac delta function and δ[·] is the Kronecker delta function. With indicator function α k = 0, H k and C k are both zeros. With α k = 1, the elements of H k and C k are independent and identically distributed (i.i.d.) Gaussian with variances ϑ H and ϑ C , respectively. We further assume that each device accesses the BS in an i.i.d. manner. Then the indicator function α k is drawn from the Bernoulli distribution as where λ is the device activity rate.
Consider an estimator to minimize the mean-square error (MSE) of H and C. It is known that the estimator which minimizes the MSE is the posterior expectation with respect to the posterior distribution [28]. Define h m ∈ C QK and c m ∈ C QK as the m-th column of H and C, respectively. Define y m ∈ C T N as the m-th column of Y. Then the posterior distribution p(H, C, α|Y) is described as are the likelihood and the prior, respectively; vector α = [α 1 , ..., α k ] T . In mMTC with a large device number K, it is computationally intractable to obtain the minimum meansquare error (MMSE) estimator. In the following section, we propose a low-complexity algorithm termed turbo message passing (Turbo-MP) to obtain an approximate MMSE solution.

IV. TURBO MESSAGE PASSING A. Algorithm Framework
The factor graph representation of p(H, C, α|Y) is shown in Fig. 3, based on which Turbo-MP is established. In the factor graph, the likelihood and prior as in (16)  treated as the variable nodes (blank circles). In specific, Turbo-MP consists of four modules referred to as Module A, B, C, and D, respectively. Module A is to obtain the estimates of the channel mean matrix H and the channel compensation matrix C by exploiting the likelihood m p(y m |h m , c m ). Modules B and C are to obtain the estimates of H and C by exploiting the priors k p(H k |α k ) and p(C k |α k ), respectively. Module D is to obtain the estimate of α k by exploiting the prior p(α k ). The estimates are passed between modules in each Turbo-MP iteration like turbo decoding [29]. For example, the output of Module B is used as the input of Module A, and vice versa. The four modules are executed iteratively until convergence. The derivation of Turbo-MP algorithm follows the sum-product rule [22] and the principle of Turbo-CS [12], [20]. In the above, y m is the subscript is the shorthand for factor node p(y m |h m , c m ); /h k,q,m denotes the set includes all elements in h m except h k,q,m . Then we define the belief of h m at factor node p(y m |h m , c m ) as Combing (17) and (18) where Σ m is the covariance matrix given by To reduce the computational complexity of channel inverse Σ −1 m , we require that A is partial orthogonal [12], i.e., Then the variance v A,post hm is given by where Σ m,i,i is the (i, i)-th element of Σ m . We note that (19) and (22)  From the sum-produce rule, the extrinsic message from factor node p(y m |h m c m ) to variable node h m is calculated as Given the Gaussian messages M(h m ) and M hm→ym (h m ), the extrinsic message is also Gaussian as where the variance v A,ext hm and the mean h A,ext m are respectively given by v A,ext The calculation to obtain the belief of c m is similar. We have the Gaussian belief M(c m ) with its mean and variance respectively given by For description convenience, the factor node p(H k |α k ) is replaced by f B k . Then the message from variable node H k to factor node p(H k |α k ) is expressed as where H B,pri Combing the Bernoulli Gaussian prior p(H k |α k ), (53), the belief of H k at factor node f B k is Bernoulli Gaussian and expressed as where λ B,post We note that the mean of vec( Define l = (m−1)Q+q. The variance of the (q, m)-th element Recall that the relationship between H k and h m is described as The above Gaussian belief assumption is widely used in message passing based iterative algorithms such as Turbo-CS [12], approximate message passing (AMP) [13] and expectation propagation (EP) [30]. Such treatment may loses some information but facilities the message updates. (Strictly speaking, the belief of the variable corresponds to a factor node instead of a Module. However, the belief of h m is associated with K factor nodes p(H k |α k ), ∀k, making the expression cumbersome. For convenience, we use the definition M B (h m ).) From the sum-product rule, we have M hm→B (h m ) = M ym→hm (h m ). Then the extrinsic message is calculated as Given the Gaussian messages M B (h m ) and M hm→B(hm) , the extrinsic message is also Gaussian with its variance and mean respectively given by Similarly to the process in Module B, each C k , ∀k in module C is estimated individually by exploiting the prior p(C k |α k ) and the messages from Modules A and D. For description convenience, f C k is used to replace factor node p(C k |α k ). The message from variable node C k to factor node p(C k |α k ) is expressed as where C C,pri k ∈ C Q×M with its (q, m)-th element c C,pri k,q,m = c A,ext k,q,m and (49), the belief of C k at factor node f C k is Bernoulli Gaussian and expressed as Then the mean of C k with respect to The variance of the (q, m)-th element of C k is given by Then we calculate the Gaussian extrinsic message with its variance and mean as follows : In Module D, we calculate the messages M α k →f B k (α k ) and M α k →f C k (α k ) as the inputs of Modules B and C, respectively. Furthermore, the device activity is detected by combing the messages from Modules B and C. According to the sumproduct rule, the message from factor node f B k to variable node α k is where Then the message from variable node α k to factor node f C k is .

(50)
Similarly to (47), the message from factor node f C k to variable node α k is where Then the message from variable node α k to factor node f B k is Define the belief of α k at factor node p(α k ) as M k (α k ). From the sum-product rule, M k (α k ) is given by .
(56) Clearly, λ D,post k (0 ≤ λ D,post k ≤ 1) indicates the probability that the k-th device is active. Therefore, we adopt a thresholdbased strategy for device activity detection aŝ where λ thr is a predetermined threshold.

F. Pilot Design and Complexity Analysis
The overall algorithm is summarized in Algorithm 1. In each iteration, the channel mean matrix H is first updated (step [2][3][4][5][6][7][8] and then the channel compensation matrix C is updated (step [9][10][11][12][13][14][15]. This is because the power of the channel mean matrix H is dominant and the iteration process effectively suppresses error propagation. Note that the prior parameters learning is shown in the next section. As mentioned in Section IV-B, the pilot matrix A is required to be partial orthogonal. To fulfill this requirement, the pilot symbols {a (t) k,n } K k=1 transmitted on the n-th subcarrier at the t-th OFDM symbol has the following property: where u i is a row vector randomly selected from an orthogonal matrix U ∈ C K×K and the selected row is different for different n, t. Combing A = [Λ 1 E 1 , ..., Λ K E 1 ] and the definition of Λ k in (4), it is easy to verify the partial orthogonality of A. We next show that when U is the discrete Fourier transform (DFT) matrix, the algorithm complexity can be further reduced. To enable the fast Fourier transform (FFT), the pilot matrix A is expressed as where S q ∈ R T N/Q×K is a row selection matrix consisting of T N/Q randomly selected rows from the K × K identity matrix. (The selected rows are different for different S q .) P = [P 1 , ..., P K ] is a column exchange matrix. In specific, in the q-th column of P k ∈ R KQ×Q , ∀k, only the k + K(q − 1)-th row is one while the others are zeros. Note that |a (t) k,n | 2 = P and matrix B has the same fast transform as A due to the relationship B = DA.
By using the FFT, the estimates of h m ∈ C QK and c m ∈ C QK , ∀m in Module A has the complexity O(QM Klog 2 K). In modules B, C and D, the estimates of H k ∈ C Q×M , C k ∈ C Q×M , and α k , ∀k, involve vector multiplications with the complexity O(QM K). As a result, the overall complexity of Turbo-MP is O(QM Klog 2 K) + O(QM K) per iteration. It is worth noting that the algorithm complexity is linear to the antenna number M and is approximately linear to the device number K.

V. PARAMETERS LEARNING
The prior parameters θ = {ϑ H , ϑ C , σ 2 w , λ} used in Turbo-MP are unknown and required to be estimated in practice. In the following, we utilize two machine learning methods, i.e., the EM and the NN approach, to learn these prior parameters.

A. EM Approach
We first use the EM algorithm [23] to learn θ. The EM process is described as where θ (i) is the estimate of θ at the i-th EM iteration. E[·|Y; θ (i) ] represents the expectation over the posterior distribution p(H, C, α|Y; θ (i) ). Note that it is difficult to obtain Similarly, the EM estimate of ϑ C is given by The EM estimate of σ 2 w is given by The EM estimate of λ is given by Turbo-MP algorithm with EM approach to learn θ is shown in Algorithm 1, which we refer to as Turbo-MP-EM. In practice, we can update H several times and then update C once to improve the algorithm stability. Besides, , , it is known the estimation accuracy of ϑ H , ϑ C , and λ by EM relies on the accuracy of the posterior approximation in Turbo-MP, and inaccurate estimation may affect the algorithm convergence. Therefore, we recommend to update ϑ H , ϑ C , and λ once after H and C is updated several times. In the first algorithm iteration, we initialize ϑ is negligible in most cases while in worst case, it continues to rise with the increase of Turbo-MP iterations. Therefore, we delete the second part of (63) in simulation.

B. NN Approach
In deep learning [3], training data can be used to train the parameters of a deep neural network. Inspired by the idea in [31] and [32], we unfold the iterations of Turbo-MP algorithm and regard it as a feed-forward NN termed Turbo-MP-NN. In specific, each iteration represents one layer of the feed-forward NN where θ is seen as the network parameter. We hope that with an appropriately defined loss function, Turbo-MP-NN can adaptively learn θ from the training data.
The structure of Turbo-MP-NN is shown in Fig. 4 which consists of L layers. Each layer has the same structure, i.e., the same linear and non-linear operations following step 2-15 in algorithm 1. To distinguish the estimates at different layer, denote H A,pri ,C A,pri and π C k obtained at the i − 1-th layer (iteration) by H A,pri,(i) , C A,pri,(i) , and π C,(i) k , respectively. In the i-th layer, the inputs consist of training data Y and the outputs of the i − 1-th layer including H A,pri,(i) , C A,pri,(i) , and {π C,(i) k } K k=1 . The loss function is defined as the normalized mean square error (NMSE) of channel estimation given by where the estimatesĤ To avoid over-fitting, we train the neural network through the layer-by-layer method [32]. Specifically, we begins from the training of the first layer, then first two layers, and finally L layers. θ is initialized following (63). For the first i layers i = 1, 2, . . . , L, we optimize θ by using the back-propagation to minimize the loss function The detailed training process is shown in Algorithm 2. Once trained, the iteration process of Turbo-MP-NN is illustrated in Algorithm 1. From the simulation, we find that for a fixed sub-block number Q, the change of the learned ϑ H and ϑ C is linear with respect to the signal-to-ratio SNR = P/σ 2 N , and the change of the σ 2 w minus σ 2 N is also linear with respect to the SNR. As such, there is no need to train the prior parameters offline at different SNR. With the loss function f (i) (θ), use back-propagation to update θ. 5: end for 6: return θ.

VI. NUMERICAL RESULTS
In this section, we evaluate the performance of the proposed Turbo-MP algorithm in channel estimation and device activity detection. The simulation setup is as follows. The BS is equipped with M = 4 or 8 antennas. P = 72 OFDM subcarriers are allocated for the random access with subcarrier spacing ∆f = 15 kHz. K = 1000 devices access the BS and transmit their signals with probability λ = 0.05 in each timeslot. We adopt the TDL-C channel model with 300 ns r.m.s delay spread and its detailed PDP can be found in TR 38.901 R14 [27]. For Turbo-MP-NN, we randomly generate 10 4 channel realizations for training and the training data is divided into minibatches of size 4. To evaluate the performance of the proposed algorithm, we use the NMSE and the detection error probability Pe = P miss + P false as the performance metrics, where P miss = 1/K k p(α k = 0|α k = 1) is the probability of miss detection and P false = 1/K k p(α k = 1|α k = 0) is the probability of false alarm. Note that for the figures showing CE performance, each data point is averaged over 5000 realizations. For the figures showing ADD performance, at each data point, the accumulative number of the detection errors is larger than 10 3 .
The baseline algorithms are as follows: • FD-GMMV-AMP: The GMMV-AMP algorithm was proposed in [18] to joint estimate the channel response on every subcarrier and detect the active devices based on the frequency-domain system model (5). We assume that the prior parameters including the channel variance on every subcarrier, noise variance σ 2 N and access probability λ are known for FD-GMMV-AMP. (The following algorithms also use σ 2 N and λ as known parameters.) • TD-Gturbo-MMV: We adopt the Gturbo-MMV algorithm [11] to achieve the ADD and CE based on the time-domain system model (7), which we refer to as TD-Gturbo-MMV. Note that the denoiser in [11] is extended and applied toH k . We further assume that the BS knows the delay spread of the time-domain channel, and thus the time-domain delay taps can be truncated to reduce the number of channel coefficients to be estimated. Specifically, there are 8 delay taps left with 4 taps at the head and 4 taps at the tail which contains 99% energy of the time-domain channel. As a Bayesian algorithm, TD-Gturbo-MMV requires the PDP of the time-domain channel as prior information. One option is to assume that the exact PDP is known. However, in practice, the exact PDP of each device is difficult to acquire. Therefore, we also consider using a flat PDP with 8 equal-power delay taps as the prior information. • TD-Vector AMP: The Vector AMP [9] algorithm is adopted as another baseline in the time-domain system model (7). Similarly, we extend and apply the denoiser in [9] toH k . TD-Vector AMP with the exact PDP and the flat PDP are both considered.   5 shows the channel estimation NMSE versus the iteration number. Both Turbo-MP-EM and Turbo-MP-NN converge and significantly outperform the other algorithms by more than 6 dB in CE NMSE. Turbo-MP-NN converges faster than Turbo-MP-EM, which implying that the neural network approach can obtain more accurate prior parameters. In addition, it is seen that FD-GMMV-AMP has a quite poor performance since the average number of the unknown variables on every subcarrier is much larger than the number of the measurements, i.e., λK T . Concerning TD-GTurbo-MMV and TD-Vector AMP, there is a performance loss when the PDP is not exactly known.  To further demonstrate the performance superiority of Turbo-MP, Fig. 6 shows the channel estimation NMSE against the SNR. With the increase of the SNR, the performance gap between Turbo-MP and the baselines becomes larger, which suggests that if the BS adopts the Turbo-MP algorithm to reach a high CE performance, the devices will consume much lower power. It is also found that the CE performance of each algorithm at different BS antenna number is similar. Fig. 7 shows the detection error probability versus the SNR. It is seen that as the antenna number M increases, the detection performances of the tested algorithms improve more than one order of magnitude. This is because the increase of BS antennas leads to a larger dimension of the block-sparsity vector. Among the tested algorithms, Turbo-MP-NN has superiority over the other algorithms especially when antenna number M = 8.    5) for Turbo-MP is only half of that (T = 10) for TD-Vector AMP and TD-Gturbo-MMV with exact PDP. It is shown that our proposed scheme can support mMTC with a dramatically reduced overhead. In Fig. 9, the detection error probability versus OFDM symbols is shown. The trend is similar to Fig. 7. Turbo-MP achieves the best ADD performance at the different number of OFDM symbols.     10 illustrates the impact of the sub-block number on the CE performance at different SNR, which implies the tradeoff between the estimation accuracy and the model accuracy. In specific, it is seen that at low SNR, a smaller sub-block number corresponds to better NMSE performance while the case is contrary at high SNR. The reason is that the NMSE performance is mainly affected by AWGN at low SNR, in which case a small sub-block number helps to improve the estimation accuracy. However, at high SNR, the CE performance is limited by the model mismatch which can be reduced by increasing the sub-block number. In practice, a proper subblock number needs to be chosen according to the wireless environment. Fig. 11 shows the miss detection probability versus false alarm probability, where we modify the threshold λ thr to reach the trade-off between miss detection and false alarm. Clearly, Turbo-MP-NN outperforms Turbo-MP-EM at different settings of sub-block number and antennas number. Besides, we see that Turbo-MP with the smaller sub-block number corresponds to a better ADD performance at SNR = −15 dB. This result is consistent with the CE performance in Fig. 10. From Fig. 10 and Fig.11, we find that the algorithms can achieve excellent device detection performance at relative low SNR while the accurate channel estimation requires a higher SNR. Such observation suggests that when the BS is equipped with multiple antennas, the bottleneck in mMTC is CE instead of ADD.
VII. CONCLUSION In this paper, a frequency-domain block-wise linear channel model was established in the MIMO-OFDM-based grantfree NOMA system to effectively compensate the channel frequency-selectivity and reduce the number with a small number of variables to be determined in channel estimation. From the perspective of Bayesian inference, we designed the low-complexity Turbo-MP algorithm to solve the ADD and CE problem, where machine learning was incorporated to learn the prior parameters. We numerically show that Turbo-MP designed for the proposed block-wise linear model significantly outperforms the state-of-the-art counterpart algorithms.