Near-Optimal Low-Complexity Hybrid Precoding for THz Massive MIMO Systems

Abstract

Terahertz (THz) communication is becoming a key technology for future 6G wireless networks because of its ultra-wide band. However, the implementation of THz communication systems confronts formidable challenges, notably beam splitting effects and high computational complexity associated with them. Our primary objective is to design a hybrid precoder that minimizes the Euclidean distance from the fully digital precoder. The analog precoding part adopts the delay-phase alternating minimization (DP-AltMin) algorithm, which divides the analog precoder into phase shifters and time delayers. This effectively addresses the beam splitting effects within THz communication by incorporating time delays. The traditional digital precoding solution, however, needs matrix inversion in THz massive multiple-input multiple-output (MIMO) communication systems, resulting in significant computational complexity and complicating the design of the analog precoder. To address this issue, we exploit the characteristics of THz massive MIMO communication systems and construct the digital precoder as a product of scale factors and semi-unitary matrices. We utilize Schatten norm and Hölder's inequality to create semi-unitary matrices after initializing the scale factors depending on the power allocation. Finally, the analog precoder and digital precoder are alternately optimized to obtain the ultimate hybrid precoding scheme. Extensive numerical simulations have demonstrated that our proposed algorithm outperforms existing methods in mitigating the beam splitting issue, improving system performance, and exhibiting lower complexity. Furthermore, our approach exhibits a more favorable alignment with practical application requirements, underlying its practicality and efficiency.

1. Introduction

Terahertz (THz) technology has emerged as a potential technique for next-generation wireless communication systems, with the goal of supporting Tbps data rates and catering to emerging ultra-high-speed applications. THz technology is well-suited for achieving high data rates due to its utilization of abundant spectrum resources [1]. THz communication, on the other hand, has obstacles such as path loss and beam splitting, making it difficult to propagate high gain THz signals over long periods of time. To address the path loss issue, researchers have leveraged massive multiple-input multiple-output (MIMO) technology to extend propagation distance in THz communication and efficiently minimize the impact of severe signal attenuation [2]. Subsequently, directional beams with high array gain are generated by precoding technology to achieve full array gain.

Hybrid precoding technology is acknowledged as a pivotal component for THz wideband mobile communication [3]. Hybrid precoding technology consists of analog RF precoding and digital baseband precoding. Traditional analog precoding, implemented with a phase shifter array, is only suitable for narrowband systems with small differences in subcarrier frequencies, resulting in beam splitting and affecting array gain in THz communication systems with significant differences in subcarrier frequencies from the center carrier frequency [4]. Furthermore, using traditional digital precoding schemes in wider antenna array systems also results in extremely high system computational complexity. Therefore, it is critical to investigate effective hybrid precoding schemes that can mitigate beam splitting while also reducing complexity [5].

In previous studies, various solutions were proposed to tackle the beam splitting problem. An antenna array structure operating in an extremely wide instantaneous bandwidth was proposed in [6], the study introduced the principles of a space-time array processor based on true-time delay (TTD) and an ultra-wideband beamformer. This structure facilitates independent control of the phase and amplitude of signals on each path. In [7], a fully connected and subarray structure was presented to enable simultaneous adjustment of amplification and phase in analog precoding. An algorithm was devised to jointly optimize phase shifters (PS) and TTD values under the constraints of each TTD device, effectively handling practical delay constraints [8]. A Delay-Phase Alternating Minimization (DP-AltMin) algorithm was proposed to address the beam splitting issue caused by the angle expansion of multipath under a channel cluster model [9]. The authors in [10] suggested a technique in which all subcarriers were projected onto a central subcarrier, aiming to simplify hardware requirements and minimize power consumption. They also constructed a generalized analog precoder based on the channel covariance matrix. However, compared to beamforming techniques based on instantaneous channel state information (CSI), this approach may result in performance degradation. A sparse RF chain antenna structure based on serial equidistant delay elements was introduced to reduce the number of RF chains and delay elements, consequently lowering hardware costs [11]. In [12], a hardware structure called true-time-delayers based delay-phase precoding (TTD-DPP) was introduced, which effectively reduces hardware costs by using some delay elements for frequency-dependent phase shifting. A joint hybrid precoding scheme based on an equivalent channel was presented in [13], which improved spectral efficiency and reduced complexity in scenarios with fewer RF chains by optimizing the hybrid precoding and composite tasks, however, the influence of beam splitting was not considered. Therefore, existing solutions either have limited performance or excessive complexity.

To address the aforementioned issues, we propose a near-optimal hybrid precoding structure with low complexity. It first performs singular value decomposition on the CSI of all subcarriers to obtain the fully digital precoder. The analog precoder is then determined using the DP-AltMin algorithm to address the impact of beam splitting. Finally, the digital precoder is designed as the product of a scaling factor and a semi-unitary matrix with a semi-orthogonal structure. By alternately optimizing the scaling factor and the semi-unitary matrix, a near-optimal digital precoding matrix is obtained. Simulation results indicate that the proposed approach offers notable benefits, including reduced computational complexity and processing time, as well as improved system performance.

This paper is structured as follows. Section II introduces the THz communication system model and channel model. In Section III, we introduce the near-optimal low-complexity hybrid precoding scheme employed in this paper. Section IV provides a simulation analysis that evaluates the system's performance and computational complexity. Finally, Section V summarizes the paper.

Notation: |⋅| denotes the magnitude of a vector, (⋅)*, (⋅)^T, (⋅)^H and ||⋅||_F denotes the conjugate, transpose, conjugate transpose and Frobenius norm of a matrix, respectively. a⊗b and a⊙b denote the Kronecker product and the Hadamard product of a and b. ∠a denotes the phase of a.

2. System Model and Channel Model

In this section, we describe the system model and channel model for the THz communication system. Table 1 presents the definitions of the channel parameters.

Table 1. Definition of parameters in the channel

2.1 System Model

We consider the delay-phase precoding (DPP) architecture [12], as shown in Fig. 1. The base station (BS) is equipped with N_t transmitting antennas and N_RF RF chains to serve users, while each user is equipped with N_r receiving antennas. The system adopts MIMO-assisted orthogonal frequency-division multiplexing (OFDM) modulation with K subcarriers [14]. For each subcarrier, the BS transmits N_s data streams to the receiver, typically with N_s ≤ N_RF << N_t. We assume that the set of transmitted data is s_k for the k-th subcarrier. The received data vector y_k can be represented as follows

Fig. 1. Delay-phase precoding architecture for THz massive MIMO.

\(\begin{align}\mathbf{y}_{k}=\sqrt{P_{t}} \mathbf{H}_{k} \mathbf{V}_{R F} \mathbf{V}_{B B} \mathbf{s}_{k}+\mathbf{n}_{k}\end{align}\),       (1)

where \(\begin{align}\sqrt{P_{t}}\end{align}\) denotes the transmit power, H_k ∈ ℂ^N_rxN_t represents the frequency-domain channel. V_BB ∈ ℂ^N_RF×N_s is the digital precoder, while ||V_RFV_BB||²_F = N_s represents the power constraint relationship for the digital precoder. s_k ∈ ℂ^N_s×1, is the transmitted data, and \(\begin{align}\mathbb{E}\left(\mathbf{s s}^{\mathrm{H}}\right)=\frac{1}{N} \mathbf{I}_{N_{s}} \cdot \mathbf{n}_{k} \in \mathbb{C}^{N_{s} \times 1}\end{align}\). n_k ∈ ℂ^N_sx1 represents the additive white Gaussian noise.

The analog precoding matrix V_RF is composed of two parts: a phase-shifting matrix that is independent of frequency and a delay matrix that is dependent on the carrier frequency. V_RF can be represented as

V_RF = V_A⊙(T_k⊗e_p),       (2)

where V_A ∈ ℂ^N_t×N_RF is a phase-shifting matrix with a constraint |V_A(i, j)| = 1. The phase-shifting matrix only adjusts the phase related to the center carrier frequency f_c. e_p ∈ ℂ^P×1 is a vector with all elements equal to 1. T_k ∈ ℂ ^MxN_RF is a frequency-dependent phase-shifting matrix generated by the delay array at the frequency point f_k. To save hardware costs, a partial connectivity approach is employed; each RF chain is connected to M delay elements. The phase shifters are divided into M groups, and each delay element is connected to Q = N_t/M phase shifters. The delay array satisfies the delay constraint T_k(m, n) = e^{-j2πf_kT(m,n)}, where T(m, n) represents the actual delay of the m-th delay element on the n-th RF chain. Each delay element T(m, n) is directly connected to Q phase shifters through the Kronecker product of T_k and e_p. The phase-shifting matrix and delay matrix work together to achieve frequency-dependent phase shift, compensating for the impact of the beam splitting effect in the wideband.

2.2 Channel Model

The design of hybrid precoders requires accurate channel information, particularly in THz communication where the channel exhibits unique characteristics, including significant path loss, severe atmospheric absorption, and limited diffraction capabilities. To accurately simulate the channel properties of THz communication systems, we employ the classical ray-tracing model. This model, which is based on geometric optics theory, is a commonly used wireless signal propagation model [15]. It simulates signal propagation in the environment by considering the paths of rays and phenomena such as reflection and refraction. In the THz frequency range, this model enables accurate simulation of signal propagation paths, attenuation, and interference. By considering crucial factors such as signal path gain, arrival angle, and delay, we can construct an accurate representation of the THz communication system channel. Building upon the ray tracing model, on the k-th subcarrier, the channel matrix H_k ∈ ℂ^N_r×N_t can be represented as

H_k = H^L_k + H^N_k,       (3)

where the channel consists of both line-of-sight (LoS) and non-line-of-sight (NLoS) components. H^L_k represents the LoS component, corresponding to single-path transmission. On the other hand, H^N_k represents the NLoS components, modeled as a cluster and involving multipath transmission [16]. Their representations are

H^L_k = 𝑎₀e^{-j2πf_kτ₀}a_r(N_r, θ^r₀)a^H_t(N_t, θ^t₀)       (4)

\(\begin{align}\mathbf{H}_{k}^{N}=\sum_{l=1}^{L_{c}} \sum_{q=1}^{L_{q}} \alpha_{l, q} e^{-j 2 \pi f_{k} \tau_{l,}} \mathbf{a}_{r}\left(N_{r}, \theta_{l, q}^{r}\right) \mathbf{a}_{t}^{\mathrm{H}}\left(N_{t}, \theta_{l, q}^{t}\right)\end{align}\).       (5)

We adopted the cluster channel model in this work, primarily focusing on the analysis of (5), where α₀ and α_{l, q} represent the complex gains for the path model and the cluster model, respectively. f_c is the center carrier frequency, and the frequency of the k-th subcarrier is denoted by \(\begin{align}f_{k}=f_{c}-\frac{B}{2}+\frac{k}{K}(k=0,1, \cdots, K-1)\end{align}\). τ₀ and τ_{l, q} represents the path delays for the single-path model and the delay for the cluster model, respectively. θ^r₀, θ^t₀ refer to the receiving and transmitting angles for the single-path model, while , θ^r_{l, q}, θ^t_{l, q} representing the receiving and transmitting angles for the cluster model. In practical channel cluster models, a clustered distribution of scattering objects causes the delays and path angles of the actual channel to have a clustered distribution; i.e., angle spread occurs [17]. An angle within the l-th cluster is modeled as θ^t_{l, q} = θ^t_l + ∆θ^t_q , where θ^t_l is the emission angle of the center cluster and ∆θ^t_q follows a Laplace distribution [18]. The vectors a_r(N_r, θ^r) and a_t(N_t, θ^t) represent the antenna array response. Using uniform linear arrays (ULAs), they can be represented as

\(\begin{align}\mathbf{a}(N, \theta)=\frac{1}{\sqrt{N}}\left[1, e^{j \pi \tau_{f_{c}}^{f_{c}} \sin \theta}, e^{j 2 \pi \frac{f_{k}}{f_{c}} \sin \theta}, \ldots, e^{j(N-1) \pi \frac{f_{k_{k}}}{f_{c}} \sin \theta}\right]^{\mathrm{T}}\end{align}\).       (6)

It can be observed that a(N, θ) is dependent on the subcarrier frequency f_k. Let us define the equivalent angle on the subcarrier as

\(\begin{align}\sin \theta_{k}=\frac{f_{k}}{f_{c}} \sin \theta, \theta \in\left[-\frac{\pi}{2}, \frac{\pi}{2}\right]\end{align}\).       (7)

2.3 Problem Modeling

Assuming a time-division duplexing (TDD) mode, the BS estimates the downlink channel through the uplink channel [19]. For a THz massive MIMO system with K subcarriers, the sum-rate R is expressed as

\(\begin{align}R\left(\mathbf{V}_{R F}, \mathbf{V}_{B B}\right)=\sum_{k=0}^{K-1} \log _{2}\left|I_{N_{r}}+\frac{P_{t}}{N_{s} \sigma^{2}} \mathbf{H}_{k} \mathbf{V}_{R F} \mathbf{V}_{B B} \mathbf{V}_{B B}^{\mathrm{H}} \mathbf{V}_{R F}^{\mathrm{H}} \mathbf{H}_{k}^{\mathrm{H}}\right|\end{align}\).       (8)

The goal is to optimize the system's sum-rate. To approximate the optimal fully digital precoder, this study discusses the search for a hybrid precoder. This involves the joint optimization of matrix variables V_RF and V_BB to address the challenge of maximizing the sum-rate. The optimization problem of precoding can be represented by

\(\begin{align}\begin{array}{l}\left(\mathbf{V}_{R F}^{o p t}, \mathbf{V}_{B B}^{o p t}\right)=\underset{\mathbf{V}_{R F}, \mathbf{V}_{B B}}{\arg \max } R\left(\mathbf{V}_{R F}, \mathbf{V}_{B B}\right) \\ \text { s.t. }\left\|\mathbf{V}_{R F} \mathbf{V}_{B B}\right\|_{F}^{2}=N_{s} .\end{array}\end{align}\)       (9)

In general, solving optimization problems with such constraints directly can be challenging. To solve the optimization problem, we can transform it into a minimum Euclidean distance problem for designing a hybrid precoder and a fully digital precoder [20]. This approximate transformation converts the original nonconvex optimization problem into a matrix factorization problem. This approach can effectively reduce the computational complexity of precoding algorithms while achieving high-gain THz wideband precoding. We convert (9) to

\(\begin{align}\min _{\mathbf{V}_{R F}, \mathbf{V}_{B B}} \sum_{k=0}^{K-1}\left\|\mathbf{F}_{o p t}-\mathbf{V}_{R F} \mathbf{V}_{B B}\right\|_{F}^{2}\end{align}\) ,       (10)

where F_opt is the optimal fully digital precoder for each subcarrier, and it can be obtained through the singular value decomposition (SVD) of the channel matrix. The singular value decomposition of H_k is given by

H_k = UAV^H,       (11)

where V and U are unitary matrices with dimensions N_t and N_r, respectively. According to the properties of SVD, the optimal fully digital precoder is formed by selecting the first N_s columns of the right unitary matrix of the channel matrix H_k, i.e.,

F_opt = [V]_{:, 1:Ns}.       (12)

To address the beam splitting effect, the analog precoding matrix adopts DP-AltMin [9]. The digital precoding matrix using traditional least square (LS) methods is extremely computationally intensive. In this paper, a semi-orthogonal structure is applied to the digital precoding matrix, representing it as a product of a scaling factor γ and a semi-unitary matrix V_DD, which are then optimized alternately. The digital precoding matrix is defined as

\(\begin{align}\mathbf{V}_{B B} \triangleq \gamma \mathbf{V}_{D D}\end{align}\),       (13)

where γ is a non-zero scaling factor, and V_DD is a semi-unitary matrix that satisfies the constraint V^H_DDV_DD = I_Ns. Thus, (13) can be represented as

V^H_BBV_BB = γ²I_Ns, γ ≠ 0,       (14)

the optimization problem in (10) can ultimately be formulated as the following constrained optimization problem

\(\begin{align}\begin{array}{l}\min _{\gamma, \mathbf{V}_{R F}, \mathbf{V}_{D D}} \sum_{k=0}^{K-1}\left\|\mathbf{F}_{o p t}-\mathbf{V}_{R F} \gamma \mathbf{V}_{D D}\right\|_{F}^{2} \\ \text { s.t. }\left\|\mathbf{V}_{R F} \mathbf{V}_{B B}\right\|_{F}^{2}=N_{s} .\end{array}\end{align}\)       (15)

3. The Proposed Hybrid Precoding Structure for THz Communication System

3.1 Analog Precoder Design

The analog precoding matrix is composed of a delay matrix and a phase shift matrix, with the following basic principle:

Given V_BB and T_k, the design of V_A according to (10) can be equivalently expressed as

\(\begin{align}\min _{\mathbf{v}_{A}} \sum_{k=0}^{K-1}\left\|\mathbf{F}_{o p t}-\left(\mathbf{V}_{\mathbf{A}} \odot\left(\mathbf{T}_{k} \otimes \mathbf{e}_{p}\right)\right) \mathbf{V}_{B B}\right\|_{F}^{2} \cdot\end{align}\).       (16)

For the objective function of (16) and utilizing the Cauchy-Schwarz inequality, it can be expressed as

\(\begin{align}\sum_{k=0}^{K-1}\left\|\left(\left(\mathbf{F}_{o p t} \mathbf{V}_{B B}^{\dagger}\right) \odot\left(\mathbf{T}_{k} \otimes e_{P}\right)^{*}-\mathbf{V}_{A}\right) \mathbf{V}_{B B}\right\|_{F}^{2}\\ \leq \sum_{k=0}^{K-1}\left\|\left(\mathbf{F}_{o p t} \mathbf{V}_{B B}^{\dagger}\right) \odot\left(\mathbf{T}_{k} \otimes e_{P}\right)^{*}-\mathbf{V}_{A}\right\|_{F}^{2} \cdot\left\|\mathbf{V}_{B B}\right\|_{F}^{2}\end{align}\).      (17)

Then, the solution of problem (16) can be determined as

\(\begin{align}\mathbf{V}_{A}=\exp \left\{j \angle\left(\sum_{k=0}^{K-1}\left\|\mathbf{V}_{B B}\right\|_{F}^{2}\left(\mathbf{F}_{o p t} \mathbf{V}_{B B}^{\dagger}\right) \odot\left(\mathbf{T}_{k} \otimes e_{P}\right)^{*}\right)\right\}\end{align}\).       (18)

Similar to the optimization process of V_A, neglecting the constant term ||V_BB||²_F, we can convert (16) into

\(\begin{align}\begin{array}{l} \min _{\mathbf{T}_{\mathrm{t}}} \sum_{k=0}^{K-1}\left\|\left(\mathbf{F}_{o p t} \mathbf{V}_{B B}^{\dagger}\right) \odot \mathbf{V}_{\mathbf{A}}^{*}-\left(\mathbf{T}_{k} \otimes \mathbf{e}_{p}\right)\right\|_{F}^{2} &=\min _{\mathbf{T}_{k}} \sum_{k=0}^{K-1}\left(\left\|\left(\mathbf{F}_{o p t} \mathbf{V}_{B B}^{\dagger}\right) \odot \mathbf{V}_{\mathbf{A}}^{*}\right\|_{F}^{2}\right. \\ & \left.-2 \operatorname{Re}\left(\operatorname{tr}\left(\left(\mathbf{F}_{o p t} \mathbf{V}_{B B}^{\dagger} \odot \mathbf{V}_{\mathbf{A}}^{*}\right)^{\mathrm{H}}\left(\mathbf{T}_{k} \otimes \mathbf{e}_{p}\right)\right)\right)+\left\|\mathbf{T}_{k} \otimes \mathbf{e}_{p}\right\|_{F}^{2}\right)\end{array}\end{align}\),       (19)

the design of T_k according to (19) can be expressed as

\(\begin{align}\max _{\mathbf{T}_{k}} \sum_{k=0}^{K-1} \operatorname{Re}\left(\operatorname{tr}\left(\left(\mathbf{F}_{o p t} \mathbf{V}_{B B}^{\dagger} \odot \mathbf{V}_{\hat{A}}^{*}\right)^{\mathrm{H}}\left(\mathbf{T}_{k} \otimes \mathbf{e}_{p}\right)\right)\right)\end{align}\),       (20)

from (20), we have D^H_k(T_k⊗e_p) = Θ^H_kT_k , where D_k = (F_optV†_BB)⊙V*_A, \(\begin{align}\Theta_{k}(m, n)=\sum_{q=1}^{Q} D_{k}((M-1) Q+q, n)\end{align}\). Substituting T_k(m, n) = e^{-j2πf_kT(m, n)} into (20), the optimization problem in (20) is found to be equivalent to the maximization problem in (21)

\(\begin{align}\max _{\mathbf{T}_{k}} \sum_{k=0}^{K-1} \operatorname{Re}\left(\operatorname{tr}\left(\boldsymbol{\Theta}_{k}^{\mathrm{H}} \mathbf{T}_{k}\right)\right)=\max _{T(m, n)} \sum_{m=1}^{M} \sum_{n=1}^{N_{k F}} \sum_{k=0}^{K-1} \operatorname{Re}\left(\Theta_{k}^{*}(m, n) e^{-j 2 \pi f_{k} T(m, n)}\right)\end{align}\).       (21)

Assuming the delay of the delay element T(m, n) is denoted as t, where t takes values within the range [0, T_max], and there are a total of S delay values. When S is sufficiently large, evenly traversing the S delay values allows us to obtain the optimal solution for t, Thus, the optimal values of T(m, n) and T_k can be obtained.

3.2 Design of the Initial Value of the Scaling Factor

Given V_RF, the goal of maximizing system performance and sum-rate, the design of V_BB can be determined as follows.

\(\begin{align}\sum_{k=0}^{K-1} \log _{2}\left|I_{N_{r}}+\frac{P_{t}}{N_{s} \sigma^{2}} \mathbf{H}_{k} \mathbf{V}_{R F} \mathbf{V}_{B B} \mathbf{V}_{B B}^{\mathrm{H}} \mathbf{V}_{R F}^{\mathrm{H}} \mathbf{H}_{k}^{\mathrm{H}}\right|\end{align}\).       (22)

According to the sum-rate formula and the water-filling power allocation principle, the optimal baseband precoding matrix is given by

V_BB = (V^H_RFV_RF)^-1/2U_eΓ_e,       (23)

where U_e represents the set of right singular vectors corresponding to the first N_s largest singular values of H_k(V^H_RFV_RF)^-1/2, and Γ_e is the diagonal matrix that represents the power allocation to each data stream using the water-filling solution. Generally, we have Γ_e ≈ I_Ns. The number of diagonal elements of V^H_RFV_RF is exactly N_t, and the off-diagonal elements can be approximated as the sum of N_t independent terms. Because of the sparsity of massive MIMO systems, the probability of having fewer independent terms than N_t is high, so it can be approximated as V^H_RFV_RF ≈ N_tI. This characteristic also proves that the optimal digital precoder for N_RF = N_s usually satisfies V^H_BBV_BB ∝ I, and the scaling factor can be further assumed to be obtained under the condition of equal power allocation for all data streams, i.e., \(\begin{align}\boldsymbol{\Gamma}_{e} \approx \sqrt{P / N_{R F}} \mathbf{I}\end{align}\). Therefore, the V_BB ≈ γU_e, and we have

\(\begin{align}\gamma=\sqrt{P /\left(N_{t} N_{R F}\right)}\end{align}\) .       (24)

3.3 Semi-Unitary Matrix and Scaling Factor Optimization for Digital Precoder Design

Given the analog precoder matrix and initial values of the scaling factors, the (15) can be further expressed as

\(\begin{align}\begin{aligned} \min _{\gamma, \mathbf{V}_{R F}, \mathbf{V}_{D D}}\left\|\mathbf{F}_{o p t}-\mathbf{V}_{R F} \gamma \mathbf{V}_{D D}\right\|_{F}^{2} & =\min _{\mathbf{V}_{p D}} \operatorname{tr}\left(\left(\mathbf{F}_{o p t}-\gamma \mathbf{V}_{R F} \mathbf{V}_{D D}\right)^{\mathrm{H}}\left(\mathbf{F}_{o p t}-\gamma \mathbf{V}_{R F} \mathbf{V}_{D D}\right)\right) \\ & =\min _{\mathbf{v}_{D D}}\left\|\mathbf{F}_{o p t}\right\|_{F}^{2}-2 \gamma \operatorname{Re}\left(\operatorname{tr}\left(\mathbf{V}_{D D} \mathbf{F}_{o p t}^{\mathrm{H}} \mathbf{V}_{R F}\right)\right)+\left\|\mathbf{V}_{R F} \mathbf{V}_{B B}\right\|_{F}^{2}\end{aligned}\end{align}\),       (25)

where the constant terms ||F_opt||²_F = N_s and ||V_RFV_BB||²_F = N_s can be obtained from the power constraint. Therefore, the above minimization problem can be rewritten as the following maximization problem

\(\begin{align}\mathbf{V}_{D D}^{o p t}=\underset{\mathbf{V}_{p D}}{\arg \max } \gamma \operatorname{Re}\left(\operatorname{tr}\left(\mathbf{V}_{D D} \mathbf{F}_{o p t}^{\mathrm{H}} \mathbf{V}_{R F}\right)\right)\end{align}\),       (26)

for (26), we have

\(\begin{align}\begin{array}{l}\gamma \operatorname{Re}\left(\operatorname{tr}\left(\mathbf{V}_{D D} \mathbf{F}_{o p t}^{\mathrm{H}} \mathbf{V}_{R F}\right)\right) & \leq \left|\operatorname{tr}\left(\gamma \mathbf{V}_{D D} \mathbf{F}_{o p t}^{\mathrm{H}} \mathbf{V}_{R F}\right)\right| \\ & \stackrel{(a)}{\leq}\left(t r\left|\mathbf{V}_{D D}\right|^{P}\right)^{1 / p}\left(t r\left|\gamma \mathbf{F}_{o p t}^{\mathrm{H}} \mathbf{V}_{R F}\right|^{q}\right)^{1 / q} \\ & \stackrel{(b)}{=}\left\|\mathbf{V}_{D D}^{\mathrm{H}}\right\|_{\infty} \times\left\|\gamma \mathbf{F}_{o p t}^{\mathrm{H}} \mathbf{V}_{R F}\right\|_{1} \\ & =\left\|\gamma \mathbf{F}_{o p t}^{\mathrm{H}} \mathbf{V}_{R F}\right\|_{1}=\sum_{i=1}^{N_{s}} \sigma_{i}, \\\end{array}\end{align}\),       (27)

where (𝑎) follows the Hölder inequality, and p > 0, 1/p+1/q = 1, while ||⋅||_∞ and ||⋅||₁ in (b) represent the L_∞ and L₁ of the Schatten norm [21], respectively. σ_i is the i-th singular value of γF^H_optV_RF, where i = 1, 2, ..., N_s.

The truncated SVD of γF^H_optV_RF is given by

γF^H_optV_RF = U₁Λ₁V^H₁,       (28)

and when (𝑎) holds, we obtain

V_DD = V₁U^H₁.       (29)

From (29), we can obtain the digital precoding matrix V_BB

V_BB = γV_DD =γV₁U^H₁.       (30)

This shows that the solution to the digital precoding matrix avoids the problem of high-dimensional matrix inversion and only involves singular value decomposition operations, resulting in a significant reduction in system computational complexity. However, the scaling factor is still not optimal and needs further optimization.

The initial value of the scaling factor can be determined according to the water-filling power allocation rule. However, due to the certain error between the optimized semi-unitary matrix V_DD and the actual γ, further optimization of γ is needed. The optimization problem in (15) is simplified to

\(\begin{align}\begin{aligned} \min _{\gamma}\left\|\mathbf{F}_{o p t}-\mathbf{V}_{R F} \gamma \mathbf{V}_{D D}\right\|_{F}^{2} & =\min _{\gamma} \operatorname{tr}\left(\left(\mathbf{F}_{o p t}-\gamma \mathbf{V}_{R F} \mathbf{V}_{D D}\right)^{\mathrm{H}}\left(\mathbf{F}_{o p t}-\gamma \mathbf{V}_{R F} \mathbf{V}_{D D}\right)\right) \\ & =\min _{\gamma}\left\|\mathbf{F}_{o p t}\right\|_{F}^{2}-2 \gamma \operatorname{Re}\left(\operatorname{tr}\left(\mathbf{V}_{D D} \mathbf{F}_{o p t}^{\mathrm{H}} \mathbf{V}_{R F}\right)\right)+\gamma^{2}\left\|\mathbf{V}_{R F} \mathbf{V}_{D D}\right\|_{F}^{2}\end{aligned}\end{align}\),       (31)

and the optimal value of γ is γ_opt from (32)

\(\begin{align}\gamma_{\text {opt }}=\frac{\operatorname{Re}\left(\operatorname{tr}\left(\mathbf{V}_{D D} \mathbf{F}_{o p t}^{\mathrm{H}} \mathbf{V}_{R F}\right)\right)}{\left\|\mathbf{V}_{R F} \mathbf{V}_{D D}\right\|_{F}^{2}}\end{align}\).       (32)

At this point, the minimum value of (31) is achieved.

The overall framework of the proposed algorithm is as follows:

Algorithm 1: Low-Complexity Optimization Algorithm

4. Simulation results

In this section, we evaluate the performance of the proposed method in THz massive MIMO systems through numerical simulation. We compare it with several other schemes, including fully digital precoding, traditional spatially sparse precoding [20], DP-AltMin [9], TTD-DPP[12], and an improved OMP algorithm [22]. The parameter settings in this paper are shown in Table 2.

Table 2. Simulation parameters

4.1 Achievable Sum-Rate Comparison

In Fig. 2 and Fig. 3, we present the relationship between system throughput and signal-to-noise ratio (SNR) for cases where the N_RF = 4 and N_RF = 8, respectively. The configuration assumes that the N_s per subcarrier is N_s = 4. Each RF chain has M = 16 delay elements connected to it, with each delay element connected to N_t/M phase shifters. The SNR is defined as SNR = 10 log₁₀ P_t/σ², where P_t is the system transmission power and σ² is the noise power. It can be seen from Fig. 2, when N_RF = 4, as SNR increases, the achievable system sum-rate of all algorithms exhibits a linear upward trend. The hybrid precoding scheme proposed in this paper achieves higher system rates at different SNRs, particularly at high SNRs, where the sum-rate achieved by the proposed scheme is closest to that of the fully digital precoder. Fig. 3 illustrates the relationship between system sum-rate and SNR when N_RF = 8. With an increase in the number of RF chains, the overall system sum-rate improves. The algorithm proposed in this paper enables a closer approximation to fully digital precoders. This is because the proposed scheme combines a delayed-phase analog precoding approach with semi-unitary matrix digital precoding algorithm, which effectively solves the beam splitting problem and improves the system sum-rate while reducing the computational complexity and enhancing the system performance.

Fig. 2. Comparisons of achievable sum-rates for different architectures with N_RF = 4.

Fig. 3. Comparisons of achievable sum-rates for different architectures with N_RF = 8.

Fig. 4 and Fig. 5 illustrate the correlation between system sum-rate and the quantity of delay elements per RF chain at low and high SNR, respectively. We gradually increased the number of delay elements per RF chain from 1 to 32. At SNR = -10 dB, Fig. 4 illustrates that the system sum-rate shows an increasing trend as the number of delay elements per RF chain increases, the proposed approach in this paper is closest to a fully digital precoder. Fig. 5 illustrates that when the SNR reaches 10 dB, the proposed algorithm achieves the greatest system sum-rate. Considering that a higher number of delay elements can lead to increased system power consumption, a balance between achievable rates and hardware cost is crucial. Therefore, we generally set the number of delay elements per RF chain to 16. Through simulation comparisons, it is evident that the proposed scheme achieves a system sum-rate that roughly approximates that of a fully digital precoder, whereas the performance of existing hybrid precoding schemes is significantly lower

Fig. 4. Achievable sum-rate versus the number of M with SNR = -10 dB.

Fig. 5. Achievable sum-rate versus the number of M with SNR = 10 dB.

In Fig. 6 we present the sum-rate against N_t. Our observations indicate that the proposed methodology exhibits superior performance compared to other hybrid precoding schemes across various values of N_t. With the increasing of N_t, all schemes show an ascending trend in achievable system sum-rate. Notably, as the number of antennas increases, the performance advantages of the proposed method become increasingly apparent due to the consideration of eliminating the impact of beam splitting. Consequently, for the prevalent deployment of massive MIMO systems, the methodology in this paper presents heightened practicality and efficacy.

Fig. 6. Achievable sum-rate versus the number of N_t.

4.2 Energy Efficiency Comparison

Fig. 7 and Fig. 8 illustrate the relationship between system energy efficiency (EE) and the number of delay elements per RF chain at low and high SNR. The SNR values are set to -10 dB and 10 dB, respectively. EE is defined as η = R/P_total, where P_total is the power consumption. P_total = P_t + P_c + N_RFP_RF + N_RFN_tP_PS + N_RFMP_TD+ P_BB, and P_t is the system calculation cost, comprising two parts: p_c = 14.1 mW/MOps, the power consumption of the digital signal processor (DSP) for every 10 million operations (MOps), and the computational complexity N_c [23]. The power consumptions of each RF chain, phase shifter, delay element, and baseband processor are P_RF = 230 mW, P_PS = 10 mW, P_TD = 100 mW, and P_BB = 300 mW, respectively. By gradually increasing the number of delay elements per RF chain from 1 to 32, we computed the system EE under various precoding schemes. The results indicate a gradual increase in system EE with the growing number of delay elements. The system EE reaches its peak when the number of delay elements is 8 and then experiences a decline. The comparison conducted indicates that the proposed approach achieves superior EE both in high and low SNR scenarios. Combining this observation with the analysis presented in Fig. 4 and Fig. 5, it becomes evident that a higher number of delay elements per RF chain leads to a higher achievable system sum-rate, a peak in system EE occurs when the number of delay elements is approximately 8.

Fig. 7. Comparisons of EE with different M at SNR = -10 dB.

Fig. 8. Comparisons of EE with different M as SNR = 10 dB.

4.3 Computational Complexity Comparison

In the low-complexity hybrid precoding scheme, assuming that the number of iterations is I_niter, the computational complexity can be expressed as follows

𝒪(KN_tN²_r + I_niterKN_sN_RFN_t + I_niterN_RFKMS + KN_sN_tN_RF), (33)

where 𝒪(KN_tN²_r) is the computational complexity of calculating F_opt for K subcarriers, 𝒪(I_niterN_RFK(N_sN_t + MS)) is the computational complexity of updating V_RF. 𝒪(KN_sN_tN_RF) is the computational complexity of the digital precoder, which mainly involves singular value decomposition of F^H_optV_RF.

To ensure an accurate comparison of the computational complexity among these schemes, we conducted simulations using identical computer hardware and software setups and measured the running time of each simulation.

Fig. 9 illustrates the system runtime, where different numbers of BS transmitting antennas are considered. To enhance the reliability of the simulation results, we conducted 50 Monte Carlo simulations to compare the runtime of the proposed scheme with that of DP-AltMin. As shown in Fig. 9, the runtimes of both algorithms increase as the number of antennas grows. However, the proposed algorithm demonstrates a lower increase rate than the other precoding scheme. This finding is particularly relevant in the context of THz massive MIMO systems that employ a large number of antennas, as the proposed algorithm effectively reduces system complexity and shortens runtime.

Fig. 9. Runtime versus N_t.

5. Conclusion

This paper introduces a new and efficient hybrid precoding scheme for terahertz massive MIMO communication systems. This method utilizes delay-phase alternating minimization hybrid precoding to address the issue of beam splitting in THz systems. Instead of using the original LS algorithm, we employ matrix singular value decomposition to obtain the fully digital precoder. Subsequently, the digital precoder is optimized to minimize the Euclidean distance between the hybrid precoder and the fully digital precoder. The Simulation results illustrate that the proposed scheme achieves superior performance and runtime efficiency compared to other THz hybrid precoding strategies. Future research can explore the extension of this method to multi-user scenarios.