Semi-Blind Receivers for Two-Hop MIMO Relay Systems with a Combined TSTF-MSMKron Coding

Due to the increase in the number of mobile stations in recent years, cooperative relaying systems have emerged as a promising technique for improving the quality of fifth-generation (5G) wireless networks with an extension of the coverage area. In this paper, we propose a two-hop orthogonal frequency division multiplexing and code-division multiple-access (OFDM-CDMA) multiple-input multiple-output (MIMO) relay system, which combines, both at the source and relay nodes, a tensor space–time–frequency (TSTF) coding with a multiple symbol matrices Kronecker product (MSMKron), called TSTF-MSMKron coding, aiming to increase the diversity gain. It is first established that the signals received at the relay and the destination satisfy generalized Tucker models whose core tensors are the coding tensors. Assuming the coding tensors are known at both nodes, tensor models are exploited to derive two semi-blind receivers, composed of two steps, to jointly estimate symbol matrices and individual channels. Necessary conditions for parameter identifiability with each receiver are established. Extensive Monte Carlo simulation results are provided to show the impact of design parameters on the symbol error rate (SER) performance, using the zero-forcing (ZF) receiver. Next, Monte Carlo simulations illustrate the effectiveness of the proposed TSTF-MSMKron coding and semi-blind receivers, highlighting the benefit of exploiting the new coding to increase the diversity gain.


Introduction
In recent years, wireless communication systems have experienced great growth in the number of users and new applications such as autonomous vehicles, smart homes, Internet of Things (IoT) and virtual/augmented reality. Compared to fourth-generation (4G) wireless systems, 5G ones offer advantages in terms of data rate, reliability, latency, energy efficiency, and mobility. To fulfill these objectives, 5G needs to operate at high frequency bands, with more base stations in a smaller area, to provide a better reliability and transmission quality to the users [1][2][3].
That explains why in the last few years, cooperative multiple-input multiple-output (MIMO) systems have attracted a lot of attention for 5G mobile networks to increase the transmission coverage area, data rates and performance of wireless communications [4]. Cooperative MIMO systems provide spatial diversity and spatial multiplexing due to the use of multiple antennas to transmit and receive signals at each node of the systems. However, individual channel estimation in a cooperative MIMO system is a fundamental Table 1. Tensor-based MIMO cooperative systems. [24] two-hop KRST PARAFAC/ PARATUCK ALS [18] two-hop KRST nested PARAFAC ALS [19] two-hop KRST nested PARAFAC KRF [16] two-hop TST nested TD ALS-KronF [26] two-hop MKRST/MKronST nested PARAFAC KRF/KronF [28] two-hop TST coupled nested TD KronF [29] three-hop TST + PARAFAC nested TD coupled SVD/ALS [22] three-hop KRST nested PARAFAC ALS/KRF [21] multi-hop TST high-order nested TD KronF [25] multi-hop KRST nested PARAFAC KRF New two-hop TSTF+MSMKron generalized-Tucker ALS-KronF/THOSVD

Ref. System Types Codings Tensor Models Receivers
We now briefly comment on the relay systems compared in Table 1 from a historical perspective. First, it is important to note that all these systems consider an amplify-and-forward (AF) protocol at the relays except the system in [26] for which the AF protocol is compared with the decode-and-forward (DF) and estimate-and-forward (EF) ones, showing that the use of these last two protocols allows significantly improving the SER performance at the cost of an additional computational complexity at the relay. From a coding point of view, the Khatri-Rao space-time (KRST) coding was firstly used in [18,19,24] for a two-hop system and then in [25] for a multi-hop system. In [22], KRST coding is combined with a rotation coding matrix for a three-hop system.
The tensor space-time (TST) coding initialy proposed in [27], in the context of point-topoint systems, was used for a two-hop system in [16], leading to a new tensor model called nested Tucker decomposition (TD) and then for a multi-hop system in [21]. In this last reference, a new tensor model called high-order nested Tucker decomposition (HONTD) was introduced. In [28], TST coding is used in a two-hop multi-relay system where the relays directly and sequentially communicate with the destination node. The sequential transmission from the relays to the destination leads to a new coupled nested TD model. In [29], TST coding is combined with a PARAFAC coding structure for a two half-duplex relays system. Two new codings, denoted MKRST and MKronST, were proposed in [26] for a two-hop system, leading to a nested PARAFAC model for the tensor of signals received at destination which is exploited to develop closed-form semi-blind receivers for joint symbol and channel estimation.
An important difference between the systems in Table 1 and the system presented in this paper concerns the a priori information needed to eliminate scaling ambiguities. Thus, our system only requires a priori knowledge of one entry of the symbol matrices, whereas all the systems in Table 1 also require knowledge of one entry or of one row of the channel matrices, which is a much more restrictive assumption. This paper proposes a new two-hop OFDM-CDMA MIMO relay system which combines a tensor space-time-frequency (TSTF) coding with a multiple Kronecker product of symbol matrices (MSMKron) at the source and relay nodes. This new coding scheme, called TSTF-MSMKron coding, can be viewed as a generalization of the codings proposed in [26,30], with the aim of increasing the diversity gain. It is established that the signals received at the relay and destination nodes satisfy generalized Tucker models whose core tensors are the coding tensors. Assuming the coding tensors are known at both nodes, the multilinear structure of tensor models is exploited to derive two semi-blind receivers for jointly estimating the symbol matrices and individual channels. Necessary conditions for parameter identifiability with each receiver are established. Extensive Monte Carlo simulations illustrate the effectiveness of the proposed TSTF-MSMKron coding and semiblind receivers. Note that our two-hop MIMO relay system differs mainly from the systems compared in Table 1 by the proposed TSTF-MSMMKron coding scheme which induces a greater diversity gain than the codings used by the systems referenced in Table 1. Another important difference lies in the consideration of frequency-dependent channels, i.e., three-dimensional channels. These assumptions lead to received signal tensors at the relay and the destination that satisfy generalized Tucker models whose essential uniqueness is ensured by the a priori knowledge of coding tensors. Scalar ambiguities can be eliminated assuming the knowledge of only one symbol per each symbol matrix. Exploiting the tensor models of received signals allows developing two types of semi-blind receivers for estimating the information symbols and the individual channels: one is iterative based on the Bi-ALS algorithm to estimate each individual channel and the Kronecker product of symbol matrices, combined with the KronF method to separate the symbol matrices, while the other one is closed form and based on the THOSVD algorithm [31], which allows simultaneously estimating each individual channel and symbol matrix. Note that unlike almost all relay systems existing in the literature which use the AF protocol, the proposed two-hop system uses the DF protocol at the relay, which greatly facilitates its generalization to the multi-hop case.
The main contributions of the paper can be summarized as follows: • A new two-hop OFDM-CDMA system that combines a TSTF coding with a multiple Kronecker product of symbol matrices (MSMKron) at the source and relay nodes is proposed. • It is established that the tensor of signals received at each hop satisfies a generalized Tucker model. • By exploiting the tensor model of the signals received at the relay and destination nodes, two semi-blind receivers are derived to jointly estimate the individual sourcerelay and relay-destination channels and transmitted symbols. • System model uniqueness and parameter identifiability conditions for each proposed receiver are analyzed. • The performance of the TSTF-MSMKron coding and the impact of design parameters on the symbol error rate (SER) are first evaluated using the zero-forcing (ZF) receiver, i.e., under the assumption of a perfect channel knowledge, by means of extensive Monte Carlo simulations. Then, the proposed semi-blind receivers are compared for symbol and channel estimation.
The rest of the paper is organized as follows. Section 2 presents tensor preliminaries. Section 3 first describes the system model, presenting the TSTF-MSMKron coding and the signals received at the relay and destination. These signals form two tensors that satisfy generalized Tucker decompositions. In Section 4, two semi-blind receivers are proposed to jointly estimate the symbol matrices and channels. Necessary conditions for parameter identifiability are derived for each receiver. In Section 6, extensive Monte Carlo simulation results are provided to illustrate the effectiveness of the proposed two-hop relay system. Section 7 concludes the paper.
Notation: scalars, column vectors, matrices, and tensors are denoted by lowercase, boldface lowercase, boldface uppercase and boldface calligraphic letters, e.g., x, x, X, and X , respectively. The transpose, complex conjugate, complex conjugate transpose, and Moore-Penrose pseudo-inverse of X are represented by X T , X * , X H and X † , respectively. We denote by x i,j the (i, j) element and by X i. (resp. X .j ) the ith row (resp. jth column) of X ∈ C I×J . The (i 1 ,. . . ,i N ) element of the N-order tensor X ∈ C I 1 ×...×I N will be written x i 1 ,...,i N . I R and I N,R represent the identity matrix of size R × R and the identity tensor of N-order and size R × R × . . . . × R, respectively.X denotes an estimate of X andX represents the matrixX after ambiguities suppression. X I 1 ×I 2 I 3 represents an unfolding of the third-order tensor X ∈ C I 1 ×I 2 ×I 3 of dimension I 1 × I 2 I 3 . The vec and unvec operators are defined by x I 2 I 3 I 1 = vec(X I 1 ×I 2 I 3 ) ∈ C I 2 I 3 I 1 ↔ X I 1 ×I 2 I 3 = unvec(x I 2 I 3 I 1 ). By slicing the third-order tensor X along each mode, we obtain three types of matrix slices called horizontal, lateral, and frontal slices, which are denoted, respectively, as follows: . The Kronecker, Khatri-Rao, and outer products are denoted by ⊗, , and •, respectively. The operator bdiag(.) forms a block-diagonal matrix from its matrix arguments, with bdiag(X ..k ) bdiag(X ..1 , . . . , X ..K ) ∈ C KI×K J , where X ..k ∈ C I×J is the kth frontal slice of X ∈ C I×J×K .
All acronyms used in the paper are summarized after Section 7.
The mode-n product between two tensors G ∈ C R 1 ×...×R n−1 ×R n ×R n+1 ×...×R N 1 ×I N 1 +1 ×...×I N and A ∈ C I n ×R n ×I N 1 +1 ×...×I N is denoted by G × n A, with n ∈ [1, N 1 ]. This product gives an N-order tensor X ∈ C R 1 ×...×R n−1 ×I n ×R n+1 ×...×R N 1 ×I N 1 +1 ×...×I N , which is defined as [30]: x r 1 ,...,r n−1 ,i n ,r n+1 ,...,r N 1 ,i N 1 +1 ,...,i N = R n ∑ r n =1 g r 1 ,...,r n−1 ,r n ,r n+1 ,...,r N 1 ,i N 1 +1 ,...,i N a i n ,r n ,i N 1 +1 ,...,i N . (2) The sum is over the second index of the tensor A, as for the mode-n product (1) between a tensor and a matrix. For example, consider the third-order tensors G ∈ C R 1 ×I 2 ×I 3 and A ∈ C I 1 ×R 1 ×I 3 . The mode-1 product X = G × 1 A is given by: We now introduce the notion of a generalized Tucker-(N 1 , N) model for an N-order tensor X ∈ C I 1 ×...×I N , with N 1 < N, which is defined as [32,33]: where S n is an ordered subset of the set {i N 1 +1 , . . . , i N }. This model can be written in terms of mode-n products as: where G ∈ C R 1 ×...×R N 1 ×I N 1 +1 ×...×I N is the core tensor, and A (n) ∈ C I n ×R n ×J n are tensor factors for n ∈ [1, N 1 ], where J n is a subset of {I N 1 +1 , . . . , I N }. For example, let us consider two factors, where the first one is a third-order tensor A (1) ∈ C I 1 ×R 1 ×I 3 and the second one is a matrix A (2) ∈ C I 2 ×R 2 . A generalized Tucker-(2, 4) model is given by: where G ∈ C R 1 ×R 2 ×I 3 ×I 4 . In scalar form, Equation (6) can be written as:

Presentation of the Proposed Two-Hop System
Consider a two-hop MIMO OFDM-CDMA system, as illustrated in Figure 1. This system is equipped with M S , M R and M D antennas at the source, relay and destination nodes, respectively. The source-relay (H (SR) ∈ C M R ×M S ×F ) and relay-destination (H (RD) ∈ C M D ×M R ×F ) channels are assumed to be flat Rayleigh fading, which is represented by third-order tensors whose coefficients are zero-mean circularly symmetric complex Gaussian i.i.d. (independent and identically distributed) random variables that are constant during at least P transmission blocks.
The decode-and-forward (DF) protocol is considered at the relay, and the transmission occurs in two hops. During the first one, the coded symbols are transmitted by the source to the relay via the channel H (SR) and decoded at the relay. During the second one, the estimated symbols are re-encoded and then re-transmitted by the relay to the destination via the channel H (RD) . Each symbol matrix S (l) = [s (l) n l ,r l ] ∈ C N l ×R l , with r l ∈ [1, R l ], n l ∈ [1, N l ], for l ∈ [1, L], is composed of R l data streams, each one containing N l information symbols. The transmission protocol is detailed in the next section which defines the TSTF-MSMKron coding. Then, in Sections 3.3 and 3.4, the tensors of signals received at the relay and the destination will be described, respectively.

TSTF-MSMKron Coding
In the proposed relay system, the coding at the source node is composed of two steps. During the first one, a multiple Kronecker product of L symbol matrices is calculated as: where N = ∏ L l=1 N l , and R = ∏ L l=1 R l . The scalar form of (8) is: with n = n (L) n l ,r l due to the multiple Kronecker product of the symbol matrix S (l) with the other matrices S (l ) , l = 1, · · · , L and l = l.
The transmission being composed of P time-slots means each symbol s (l) n l ,r l is repeated R l times, which implies an increase of time and code diversities when increasing the dimensions N l and R l , respectively. During the second step, the MSMKron coding is combined with a tensor space-timefrequency (TSTF) coding [30] carried out by means of the (L+3)-order tensor G (S) ∈ C M S ×R 1 ×...×R L ×F×P in such a way that the tensor of signals coded at the source satisfies an (L+3)-order Tucker model given by: Note that the core tensor of this decomposition is the coding tensor G (S) . In scalar notation, the coded signals transmitted by the m th S antenna at the source, using the f th subcarrier, during the p th time slot are given by: The TSTF-MSMKron coding increases spacetime-frequency diversity, as will be illustrated in the simulations.

Tensor of Signals Received at the Relay
In the noise-free case and assuming a flat Rayleigh fading propagation channel, the signal x (SR) m R ,n 1 ,...,n L , f ,p received at the m th R antenna of the relay, during the n th l symbol period of the p th block and associated with the f th subcarrier, is given by: In terms of mode-n products, we have: Note that the transmission via channel H (SR) can be interpreted as a mode-1 linear transformation applied to the tensor V (S) of coded signals. Substituting (11) into (12) gives the signal received at the relay written in scalar form as: The signals received at the relay form the tensor X (SR) that satisfies a generalized Tucker-(L + 1, L + 3) model given by: where S (l) represents the symbol matrices encoded by the TSTF-MSMKron coding for l ∈ [1, L], and G (S) is the core tensor of the Tucker model. As is well known, knowledge of the core tensor implies the uniqueness of this model. Combining modes 2 to L + 1 of tensors G (S) and X (SR) results in contracted forms G (S) c ∈ C M S ×R×F×P and X (SR) c ∈ C M R ×N×F×P , and Equation (15) can be rewritten as: From the Tucker model (16), it is easy to deduce the following matrix unfoldings of the tensor X (SR) : with H (SR) .. f ∈ C M R ×M S and bdiag(.) previously defined in the notation. Note that the (17) is associated with FP repetitions of the symbol matrices inducing time-frequency diversity for the system.
The block structure of the matrix unfoldings G  (17) and (19), respectively, is defined as follows: (20) is a block-diagonal matrix, formed of F diagonal blocks of dimension PR × M S , each block being formed of M S column vectors corresponding to a vectorized form of the tensor slice G To illustrate the matrix unfolding (21), consider the case where R = P = M S = F = 2, leading to the following matrix:

Tensor of Signals Received at the Destination
With the DF protocol, the symbols received at the relay are first decoded by means of one of the receivers described in Section 4, leading to the estimated symbol matriceŝ S (l) , which are also written as S (l) R . The estimated symbols are then re-encoded at the relay using a TSTF-MSMKron coding, with the tensor coding G (R) ∈ C M R ×R 1 ×...×R L ×F×P . The reencoded signals are transmitted by the relay to the destination via the channel H (RD) ∈ C M D ×M R ×F . The signals received at the destination are similar to the signals received at the relay, defined by Equations (14) and (15), with the following correspondences: Similar to (14), in the noiseless case, the signal received at the m th D antenna of the destination node, during the n th l symbol period of the p th time block and associated with the f th subcarrier, is given by: and the generalized Tucker-(L + 1, L + 3) model (15) becomes: where X (RD) ∈ C M D ×N 1 ×...×N L ×F×P . Matrix unfoldings of this tensor can be deduced from (17)- (19) using the correspondences (23) and (24)  The system design parameters and the definitions of the system matrices and tensors are summarized in Tables 2 and 3, respectively.   Table 3. System matrices and tensors.

Semi-Blind Receivers
In this section, two semi-blind receivers are proposed to estimate the channel tensors and symbol matrices at the relay and destination nodes. We assume that the coding tensors G (S) and G (R) are known. We also assume that one symbol of each symbol matrix is known to eliminate scalar ambiguities. The symbol matrices and the channel tensor H (SR) are estimated at the relay, while the symbol matrices and the channel tensor H (RD) are estimated at the destination. The proposed receivers are detailed for the relay. The same receivers can be derived for the destination, using the correspondences (23) and (24). The first one is based on the alternating least squares (ALS) algorithm to estimate the channel and the Kronecker product of symbol matrices, which is followed by the Kronecker factorization (KronF) method to separate the symbol matrices, while the second one is a closed-form solution allowing to jointly estimate the channel and the symbol matrices by means of the truncated higher-order singular value decomposition (THOSVD) algorithm.

Bi-ALS-KronF Receiver
In the first step, the bi-alternating least squares (Bi-ALS) algorithm is used to jointly estimate the MSMKron product S and the channel tensor H (SR) . Then, the KronF algorithm is applied to separate the symbol matrices. The Bi-ALS algorithm results from the minimization of the following cost function deduced from Equation (16): where · F is the Frobenius norm. The Bi-ALS method replaces the optimization problem (27) by two LS sub-problems deduced from the matrix unfoldings (17) and (18), leading to the alternate minimization of the following LS criteria: min S X (SR) The update equations at iteration [it] are given by: The matrices (I FP ⊗ S)G  (30) and (31) as: The Bi-ALS algorithms (32) and (33) are simplified versions of (30) and (31) in terms of pseudo-inverses computation at the price of additional constraints on the design parameters.
The error at the [it] th iteration, deduced from (17), is considered for deciding the convergence of the Bi-ALS algorithm: Convergence at the [it] th iteration is declared when this error does not significantly change between two successive iterations, i.e., |err[it − 1] − err[it]| ≤ , where is a predefined threshold. Since the core tensor G (S) is assumed to be known, there is no permutation ambiguity, and the generalized Tucker model (16) is unique up to scalar scaling ambiguities. The LS estimatesĤ (SR) FM S ×M R andŜ, at convergence, after correcting the ambiguities are given by: For eliminating these scaling ambiguities, it is sufficient to assume that one element of S is known a priori, e.g., s 11 = 1. Under this assumption, λ (S) is calculated as: λ (S) =ŝ 11 . The symbol matrices S (l) are then estimated by means of the KronF algorithm presented in Appendix A, minimizing the following LS cost function: After applying the KronF algorithm, the estimated symbol matrixŜ (l) is obtained by unvectorizingŝ (l) as:Ŝ (l) = unvec(ŝ (l) ) ∈ C N l ×R l , and assuming s (l) 11 = 1, the scalar ambiguity is corrected by: As mentioned previously, the Bi-ALS-KronF receiver at the destination can be deduced from the one at the relay, using the correspondences (23) and (24), to estimate the channel H (RD) ∈ C M D ×M R ×F and the symbol matrices denoted S (l) R ∈ C N l ×R l . To eliminate the scaling ambiguities in the second hop, we use the same relation (38) for the KronF algorithm. At each hop, the estimated symbols are obtained after a projection onto the symbol alphabet. The Bi-ALS-KronF algorithm is summarized in Algorithm 1.

THOSVD-Based Receiver
The THOSVD-based receiver is proposed to jointly estimate the channels and the symbol matrices. This closed-form solution can be viewed as a generalization of the KronF algorithm used to separate the symbol matrices. The difference is that we can now simultaneously estimate all the matrices (H   (8) and (21), we deduce the following LS estimate of the multiple Kronecker product: with Z (SR) ∈ C M R N×FM S R . The unfolding G are jointly estimated by means of the rank-one approximationbased KronF algorithm, which is described in Appendix A. The THOSVD receiver at the destination is deduced from the one at the relay, using the correspondences (23) and (24)

Zero-Forcing (ZF)-KronF Receiver
To evaluate the impact of the design parameters on the system performance, we use the zero-forcing (ZF)-KronF receiver, which assumes a perfect channel knowledge. The LS estimate of S is obtained using (31) As for the Bi-ALS algorithm, the use of (40) or (41) implies the following necessary conditions: R ≤ PFM R or R ≤ PFM S , and M S ≤ M R . Then, the symbol matrices S (l) are estimated using the KronF algorithm as in the second step of the Bi-ALS-KronF receiver. For the second hop, the ZF-KronF receiver is similar to the one in the first hop with the correspondences (23) and (24) Table 4 summarizes the necessary conditions for parameter identifiability with each receiver. Comparing the identifiability conditions for the Bi-ALS-KronF algorithms (32) and (33) with the ones for the Bi-ALS-KronF algorithms (30) and (31), we can deduce some implications. Indeed, for the estimate (32), the conditions M S ≤ PR and R ≤ N imply M S ≤ PN, i.e., the identifiability condition for the LS solution (30). For the estimate (33), the conditions R ≤ PFM S and M S ≤ M R imply R ≤ PFM R , i.e., the identifiability condition for the LS solution (31). In other words, if the identifiability conditions for (32) and (33) are satisfied, then the ones for the Bi-ALS algorithm (30) and (31) are automatically satisfied. Note also that R ≤ PFM S and M S ≤ PR imply R ≤ P 2 FR, which is always satisfied. Therefore, the condition M S ≤ PR can be discarded. We can also conclude that the THOSVD receiver is more restrictive than the Bi-ALS receivers in the sense that a higher value of P is required, implying a reduction in the transmission rate. As the ZF-KronF receiver (41) only esti-mates the symbol matrices, its identifiability conditions are a subset of those of the second Bi-ALS-KronF receiver. Table 4. Identifiability conditions for the receivers.

Receiver
Identifiability Conditions (First Hop)

Computational Complexity
In this section, we compare the computational complexity of the proposed THOSVD and Bi-ALS-KronF receivers by evaluating the cost of SVD calculation, which is the most expensive matrix operation. Note that for a matrix of dimensions I × J, the complexity of SVD computation is O(I J min(I, J)). The complexities are evaluated by taking the identifiability conditions into account.
The computational complexity of the HOSVD algorithm for an N-th-order tensor X ∈ R I 1 ×···×I N is of the order of O (∑ N n=1 I n ) ∏ N q=1 I q if I n ≤ ∏ N q =n I q , requiring to compute N SVDs of I n × I n+1 . . . I N I 1 . . . I n−1 matrices for n ∈ [1, N].
The ALS algorithm requires, at each iteration, the overall computational complexity O R 2 ∑ N n=1 (∏ N q =n I q ) to compute the PARAFAC decomposition of a tensor X ∈ R I 1 ×···×I N assumed to be of rank R. This algorithm requires calculating N LS estimates, which needs to pseudo-inverse ∏ N q =n I q × R matrices, for n ∈ [1, N]. For estimating the L symbol matrices from their Kronecker product, the KronF algorithm has a complexity of O((∑ L l=1 N l R l ) ∏ L q=1 N q R q ) flops. In Table 5, the computational complexities of the Bi-ALS-KronF and THOSVD receivers are compared for the first hop. The computational complexities for the second hop can be easily derived using the correspondences (24) between the dimensions. Table 5. Computational complexity of the Bi-AKS-KronF and THOSVD algorithms at the first hop.

Algorithms Computational Complexity
Bi-ALS-KronF (30) and Bi-ALS-KronF (32) and Note that simplifying the pseudo-inverses in (30) and (31) results in less computational complexity for the Bi-ALS-KronF (32) and (33) than for Bi-ALS-KronF (30) and (31). Regarding the computational complexity of the closed and form THOSVD and based receiver, it is generally lower than the one of the iterative Bi-ALS algorithms, which depends on the number of iterations needed for convergence.

Simulation Results
In this section, we evaluate the performance of the proposed two-hop OFDM-CDMA MIMO system and the associated receivers. First, in Section 6.1, we describe the simulations and present the considered performance criteria. In Section 6.2, we study the impact of design parameters on the symbol error rate (SER), using the ZF-KronF receiver. Finally, in Section 6.3, the proposed semi-blind receivers are compared in terms of SER and channel normalized mean square error (NMSE).

Description of the Simulations
The noisy signals received at each hop, Y (SR) and Y (RD) , respectively, are simulated as: where N (SR) ∈ C M R ×N 1 ×...×N L ×F×P and N (RD) ∈ C M D ×N 1 ×...×N L ×F×P are additive white Gaussian noise (AWGN) tensors whose entries are zero-mean circularly symmetric complexvalued Gaussian random variables, the tensors X (SR) and X (RD) contain the noise-free received signals obtained by means of Equations (15) and (26), respectively, and α (SR) and α (RD) allow fixing the signal-to-noise ratio (SNR) calculated as: which gives α (SR) = X (SR) , are randomly generated from the 16-QAM (Quadrature Amplitude Modulation) alphabet with a uniform distribution. It is worth mentioning that our proposed coding scheme and semi-blind receivers are not dependent on a specific choice for the modulation format as presented in [34,35]. The proposed system may operate with any modulation, although the resulting SER performance and transmission rate will be affected by this choice. For instance, increasing the modulation cardinality of M-PSK (phase-shift keying) or M-QAM type constellations (under the same total transmit power constraint) would result in a higher transmission rate at the cost of an SER performance degradation. In this work, we have adopted 16-QAM since it offers a good tradeoff between SER performance and transmission rate As mentioned before, the coding tensors are designed for each Monte Carlo run: in such a way that, their matrix unfoldings G (S) PFM R ×R are truncated DFT matrices. The performance criteria, plotted versus SNR, are calculated as: whereẐ k is the tensor Z k estimated at the k th run, with Z k ∈ {H }. The SER and NMSE are calculated by averaging the results over K = 5.10 4 Monte Carlo runs, after truncating the 5% worse and 5% better values to eliminate the influence of illconvergence and outliers.
The transmission rate T is given by: where ∑ L l=1 N l R l corresponds to the total number of transmitted symbols, L is the number of symbols assumed to be a priori known for ambiguity suppression, and µ denotes the number of constellation points. Note that increasing the number N l of symbols in the symbol matrix S (l) induces an increase of coding diversity and a lower transmission rate T, while an increase of the number P of repetitions implies a decrease of T.

Impact of Design Parameters
In this section, we evaluate the SER performance of the proposed system under perfect channel knowledge. In this case, we use the ZF-KronF receiver to estimate the transmitted symbol matrices by means of Equation (41). The results presented in Figures 4-9 were obtained for both hops, but due to lack of place, some SER results are shown only for the relay. All parameters used for the simulations are provided in Table 6. Note that the default values of these parameters are chosen equal to two. The corresponding transmission rates are given in Table 7. Figure 3 shows the impact on the SER for different numbers of symbols per data stream: N 1 = N 2 ∈ {8, 12, 16}, where S relay and S dest denote the SER at the relay and the destination, respectively. From these simulation results, it can be concluded that the SER is improved when the numbers of symbols increase, which implies an increase of coding diversity, since N = N 1 N 2 is a dimension of the contracted form Y (SR) c and Y (RD) c of the data tensors, which is not the case for R = R 1 R 2 . On the other hand, the transmission rate decreases as shown in Table 7. In addition, note that the SER at the relay is better than the one at the destination. This happens because with the DF protocol, the symbols are estimated and decoded before they are retransmitted by the relay to the destination, which induces a propagation error due to the decoding. Table 6. Parameters for the simulations.            (32) and (33) and ZF receivers at the relay.  (32) and (33) and ZF receivers at the destination.  Figure 4 compares the SER for three different data stream numbers: R 1 = R 2 ∈ {4, 6, 8}. From this figure, it can be concluded that increasing R 1 and R 2 implies an increase of the number of symbols to be estimated without increasing the number of data in the tensor Y (SR) for performing the symbol estimation, thus inducing a degradation of the SER, while the transmission rate increases (see Table 7).

Impact of Parameters
In Figure 5, the simulation results compare the SER global with the individual SERs for S (1) and S (2) when N 1 = 4, N 2 = 12 and R 1 = R 2 = 2. For this configuration, the Kronecker product between S (1) and S (2) induces a greater diversity for S (1) than for S (2) due to the fact that each symbol of S (1) is repeated 12R 2 times while each symbol of S (2) is repeated only 4R 1 times. That implies an SER smaller for S (1) than for S (2) . Figure 6 presents the results considering different configurations for the numbers of subcarriers (F) and time blocks (P). Note that a performance improvement is obtained when F and/or P are/is increased due to an increase of frequency and/or time diversities. On the other hand, the transmission rate decreases. We can also remark that for the same value of the product FP = 8 or FP = 16, the diversity gain is the same, implying very close SERs, which illustrates the symmetric role played by the frequency and time diversities in the SER performance.
In Figure 7, we compare the SER for different numbers of symbol matrices (L ∈ {2, 3, 5}). The design parameters have been chosen so that the transmission rate is the same for the three values of L. The MSMKron scheme with L = 5 provides the best SER performance in comparison with L ∈ {2, 3}. These results corroborate the coding gain provided by the Kronecker product of symbol matrices.
In Figure 8, the SERs are plotted for different configurations of antenna numbers (M S , M R , M D ) ∈ { (2,4,6), (2,4,2), (4, 2, 6)}. Comparing these configurations, we note that the best SER is obtained when M D > M R > M S . For the configuration (4, 2, 6), the SER is not good both at the relay and the destination, because the identifiability condition (M S ≤ M R ) at the relay is not satisfied. For the configuration (2,4,2), the SER at the relay is similar to the one for the configuration (2, 4, 6) because the antenna numbers (M S , M R ) are the same for both configurations, but the SER at the destination is not good because the identifiability condition (M R ≤ M D ) at the destination is not satisfied for the configuration (2, 4, 2), which is not the case of the configuration (2,4,6). With this last configuration, we note that the SER at the relay is better than the one at the destination.
In Figure 9, the proposed TSTF-MSMKron coding is compared with the TSTF coding, i.e., using a single symbol matrix S ∈ C N×R instead of a multiple Kronecker product of symbol matrices. With the TSTF coding, the symbol matrix is estimated using Equation (31), and the transmission rate is given by: For both codings, the number (14) of transmitted symbols is the same. See the design parameters in Table 6.
As expected, from Figure 9, we conclude that the TSTF-MSMKron coding gives a better SER than the TSTF coding thanks to a greater coding diversity brought by the Kronecker product of symbol matrices. As a counterpart, the transmission rate with the TSTF-MSMKron coding is smaller than the one with the TSTF coding. See Table 7.

Comparison of THOSVD and Bi-ALS-KronF Receivers
In the next experiments, we compare the SERs obtained with the proposed semiblind and ZF-KronF receivers. First, the results are presented in terms of SER at the relay (S relay - Figure 10) and the destination (S dest. - Figure 11). Then, we compare the performance of semi-blind receivers in terms of channel NMSE at each hop ( Figure 12). For these simulations, the design parameters are fixed with the following values: M S = 2, M R = M D = 4, N 1 = N 2 = 4, R 1 = R 2 = 2, P = 18, and F = 2.
From Figures 10 and 11, we can conclude that the THOSVD receiver provides a better SER performance than the Bi-ALS-KronF receiver. That is due to the closed form of the THOSVD receiver allowing to jointly estimate the channel and symbol matrices, while the Bi-ALS-KronF receiver is composed of two steps, one iterative and one closed form. On the other hand, the THOSVD receiver is more constraining in terms of identifiability conditions (M S R ≤ P) than the Bi-ALS-KronF receiver, inducing a reduction of the transmission rate, as can be seen in Table 7. It can also be noted that the SER at the relay is better than the one at destination due to the error propagation caused by decoding at the relay. As expected, the ZF-KronF receiver provides the best SER due to an a priori knowledge of the channels.
In Figure 12, the channel NMSE results obtained at each hop are plotted. Note that the THOSVD receiver gives better results than the Bi-ALS-KronF one. As for the SER, this is because the THOSVD is a closed-form solution, while the Bi-ALS algorithm is iterative. Moreover, the channel estimation in the first hop is slightly better than the one in the second hop. This is due to error propagation in the re-transmission of symbol matrices after decoding at the relay.
Note that considering non-coherent receivers [36][37][38] would imply avoiding the assumption about the knowledge of the coding tensors G (S) and G (R) used at the source and the relay, which would require a fully blind approach. Such a non-coherent assumption would destroy the essential uniqueness property of the estimated channels and symbols (up to scaling ambiguities). More specifically, in the non-coherent case, the Tucker models defined in Equations (15) and (26) would be affected by rotational ambiguities, which means that the channel tensors and symbol matrices estimated at the relay and destination nodes would be linked to the true ones via a transformation by a nonsingular matrix. It should be mentioned that one possible way to ensure the successful decoding of transmitted symbols in the non-coherent case, where such rotational ambiguities are present, is to consider Grassmannian constellations for symbol matrices, as proposed in [39,40].

Conclusions
In this paper, we have proposed a new two-hop CDMA-OFDM MIMO system which combines a tensor space-time-frequency (TSTF) coding with a multiple Kronecker product of symbol matrices, leading to the so-called TSTF-MSMKron coding. This new coding makes it possible to improve the gains in diversity and throughput. We have shown that the tensors of signals received at the relay and destination nodes satisfy two generalized Tucker models whose core tensors are the coding tensors.
Assuming these coding tensors are known, two semi-blind receivers have been derived to jointly estimate the transmitted information symbols and the channels. One, called the Bi-ALS-KronF receiver, is composed of two stages. In the first stage, the iterative ALS algorithm is used to estimate the channel and the Kronecker product of symbol matrices, while in the second stage, the KronF method is applied to separate the symbol matrices. The other one, called THOSVD receiver, is a closed-form solution which allows simultaneously estimating the channel and the symbol matrices by means of SVD computations as with the KronF method. Necessary conditions for system identifiability have been established for each receiver, showing that the THOSVD receiver is more constraining than the Bi-ALS-KronF one for the choice of the number of time blocks and consequently from the data rate point of view.
It is worth mentioning that the proposed two-hop system can be easily extended to the multi-hop case owing to the use of the DF protocol at the relay, since the tensor models for the signals received at the relays and destination have the same structure (generalized Tucker models), with the correspondences (23) and (24) established between the first and second hops. These correspondences can be easily generalized to more than two hops if the same coding scheme is used at each relay.
Extensive Monte Carlo simulations have allowed illustrating the impact of all the design parameters on the SER performance using the ZF receiver. In particular, the diversity gain brought by each parameter of the TSTF-MSMKron coding has been analyzed. The performances of the proposed semi-blind receivers have been compared in terms of SER and channel NMSE. As expected, the THOSVD closed-form receiver outperforms the iterative Bi-Als-KronF receiver. Moreover, a comparison with the standard TSTF coding has corroborated the SER improvement brought by the MSMKron coding, which allows increasing the diversity gain.
Note that we have numerically evaluated the SER performance under different schemes, assuming 16-QAM constellation for all the symbol matrices involved in our MSMKron coding scheme. At this point, we do not have a theoretical SER performance evaluation. Deriving an analytical Cramer-Rao bound (CRB) for the estimated channels and symbols is challenging, and this constitutes an important perspective for this work.
Among some other perspectives of this work, we can mention an extension of the proposed relaying system to the multi-hop case using the amplify-and-forward (AF) protocol and taking resource allocation tensors into account. Such considerations will lead to new tensor models and therefore new semi-blind receivers. Other extensions concern the development of relaying systems with TSTF-MSMKron coding for double-directional dual-polarized MIMO systems and intelligent reflecting surfaces (IRS)-assisted systems.

Appendix A. Kronecker Factorization (KronF) Algorithm
In this section, the KronF algorithm is presented for estimating the matrix factors of a multiple Kronecker product C = A (1) ⊗ . . . ⊗ A (N) ∈ C I 1 ...I N ×R 1 ...R N , with A (n) ∈ C I n ×R n , for n ∈ [1, N], by minimizing the LS cost function: min A (n) ,n∈ [1,N] Following the idea introduced in [41] for a Kronecker product of two matrices, the problem is solved by rewriting the cost function (A1) in terms of approximation of a rank-one tensor built as the outer product of vectorized forms of the matrix factors, as: min a (n) ,n∈ [1,N] C − a (1) where a (n) = vec(A (n) ) ∈ C R n I n , and C ∈ C R 1 I 1 ×...×R N I N is the rank-one tensor obtained by reshaping the multiple Kronecker product: C = reshape(C, [R 1 I 1 , . . . , R N I N ]).
Each vector a (n) is estimated using the THOSVD algorithm, and the matrix factor estimateÂ (n) is deduced using the unvec operator [29,32,42], with a scalar ambiguity which can be eliminated assuming the knowledge of one element of A (n) , e.g., a (n) 11 = 1, which leads to the following corrected estimate: (A4)