{'fraction_non_alphanumeric': 0.03298297293075361, 'fraction_numerical': 0.007185735645265671, 'mean_word_length': 4.675997466751109, 'pattern_counts': {'":': 0, '<': 2, '<?xml version=': 0, '>': 0, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 7, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Spatially Sparse Precoding in Wideband Hybrid Terahertz Massive MIMO Systems

In terahertz (THz) massive multiple-input multiple-output (MIMO) systems, the combination of huge bandwidth and massive antennas results in severe beam split, thus making the conventional phase-shifter based hybrid precoding architecture ineffective. With the incorporation of true-time-delay (TTD) lines in the hardware implementation of the analog precoders, delay-phase precoding (DPP) emerges as a promising architecture to effectively overcome beam split. However, existing DPP approaches suffer from poor performance, high complexity, and weak robustness in practical THz channels. In this paper, we propose a novel DPP approach in wideband THz massive MIMO systems. First, the optimization problem is converted into a compressive sensing (CS) form, which can be solved by the extended spatially sparse precoding (SSP) algorithm. To compensate for beam split, frequency-dependent measurement matrices are introduced, which can be approximately realized by feasible phase and delay codebooks. Then, several efficient atom selection techniques are developed to further reduce the complexity of extended SSP. In simulation, the proposed DPP approach achieves superior performance, complexity, and robustness by using it alone or in combination with existing DPP approaches.


I. INTRODUCTION
In future 6G wireless communication systems, it is foreseeable that ultra-high data rates will be required to support many cutting-edge applications, such as virtual reality, augmented reality, and digital twins [1]. With abundant bandwidth resources, terahertz (THz) is a promising frequency band to meet the sharply increased data rate requirement [2], [3]. Specifically, the bandwidth provided by the THz band could reach up to several tens of GHz, which is much larger than that in the millimeter-wave (mmWave) band used in 5G [2]. However, higher carrier frequencies make the THz signals suffer from more severe attenuation than signals in lower frequency bands [4]. To alleviate this issue, THz is usually combined with massive multiple-input multiple-output (MIMO), where directional beams with high array gains can be generated by massive antennas at transceivers [4]- [6]. Furthermore, to relieve the manufacturing cost and power consumption, it is a common practice to adopt the hybrid precoding architecture in THz massive MIMO systems, where precoding is decomposed into a high-dimensional analog precoder and a lowdimensional digital precoder [7], [8]. Thanks to the sparsity of THz channels in the angular domain, satisfactory performance can still be achieved with only a few radio-frequency (RF) chains using hybrid precoding. Next, we will briefly review prior works about hybrid precoding in THz massive MIMO systems. Since mmWave channels and THz channels have many similar properties, some works based on the mmWave model are included as well.

A. Prior Works
In most prior works, the analog precoder is implemented by phase shifters. Since the optimal solutions are hard to obtain in the presence of the non-convex unit modulus constraint of phase shifts, the approximate orthogonal properties of the optimal analog and digital precoders are exploited in [9] to simplify the optimization process. In the well-known spatially sparse precoding (SSP) algorithm [10], the sum rate maximization problem is first converted to an approximately equivalent while much simpler matrix decomposition problem, where the distance between the fully-digital precoder and the hybrid precoder is minimized. Then, by exploiting the structures of optimal precoders, the optimization problem is further converted to a typical sparse recovery problem that can be readily solved by various compressive sensing (CS) algorithms. In [11], the complexity of SSP is reduced by using orthogonal codebooks and beamspace singular value decomposition (SVD). To achieve better performance without constraining the form of the analog precoder, an alternative minimization (AM) algorithm is proposed in [12] to optimize the analog and digital precoders alternately. In [13], [14], deep neural networks are used to predict hybrid precoders based on channel state information (CSI). Without the need of estimating the high-dimensional full CSI, codebook based approaches are appealing in practice, where efficient codebooks and beam sweeping procedures are developed to reduce the overhead [15], [16].
Although the above phase-shifter-based hybrid precoding architecture works well in most narrowband systems [8], it will suffer from severe performance degradation in wideband systems, especially when the number of antennas and the fractional bandwidth (ratio between bandwidth and central frequency) are large, e.g., THz massive MIMO. The reason is the existence of the beam split effect 1 , which is proportional to both the number of antennas and the fractional bandwidth [17]- [19]. With beam split, there is a natural contradiction, such that the same physical channel path will have different equivalent angles at different subcarriers and make the optimal fully-digital precoders frequency-dependent while the analog precoder implemented by phase shifters is shared by all subcarriers and frequency-independent [19].
To minimize the impact of beam split, several methods have been proposed from the perspective of algorithm redesign [19]- [23]. In [19], a beam broadening approach has been proposed where wider beams are constructed through subarray coordination to cover the split directions at all subcarriers. Since the AM algorithm works with general channels including those with beam split, it is used in [20] to solve the wideband hybrid precoding problem. More generally, AM is used in [21] to develop mmWave and THz hybrid precoding algorithms for various analog architectures, such as the quantized phase shifters and the use of switches for dynamic connections. In [22], [23], codebooks are revised to maximize the minimum array gain and the total array gain achieved across all subcarriers, respectively. Nevertheless, these methods still use a frequency-independent analog precoder, therefore will lose most of the effectiveness when beam split is relatively severe.
Since beam split essentially comes from the delay of electromagnetic wave travelling through the antenna array [17], a more direct and thorough solution is to compensate for the delay from the perspective of hardware redesign. In the emerging delay-phase precoding (DPP) architecture, true-time-delay (TTD) lines are incorporated apart from phase shifters to implement analog precoders. Since the phase shift of a TTD line is proportional to both its delay parameter and the subcarrier frequency, frequency-dependent analog precoders can be realized. Assuming the path channel model, heuristic DPP algorithms that exploit the angles of channel paths sorted by their path gains are proposed in [19], [24], which can achieve near-optimal performance with low complexity. To support cluster channels, the AM algorithm is used in [25] to optimize the time delays, phase shifts, and digital precoders alternately. To improve energy efficiency, fixed TTD lines are used along with a switch network in [26]. In [27], multi-resolution time-delay codebooks are designed through subarray coordination for efficient beam sweeping, while the level of beam split is controlled by delay-phase configuration to accelerate beam tracking in [28]. In [29], the effectiveness of the DPP architecture is validated by practical 6G THz massive MIMO channels.
Among existing DPP algorithms, it is hard to find one that satisfies both the strict requirements for performance and complexity in 6G communication systems with practical THz channels.
Specifically, the performance of heuristic DPP algorithms will degrade in cluster channels while iterative optimization based DPP algorithms usually suffer from high complexity. The SSP algorithm originally designed for narrowband systems has many attractive characteristics thanks to the exploitation of precoder structures [10]. Nevertheless, its extension to wideband DPP with beam split considered is not investigated yet.

B. Contributions
In this paper, we first propose a novel extended SSP algorithm to solve the DPP problem in wideband THz massive MIMO systems and then propose several tailored techniques to further reduce the algorithm's complexity. The main contributions are as follows: • Following the key idea of SSP, we first formulate the wideband DPP problem as a matrix decomposition problem to approximate fully-digital precoders with hybrid precoders at all subcarriers. Then, the optimization problem is further formulated as a multiplemeasurement-vector (MMV) CS problem by exploiting the structures of optimal precoders. Eventually, the simultaneous orthogonal matching pursuit (SOMP) algorithm is used to solve the MMV CS problem.
• To recover the common sparsity structure destroyed by beam split and sustain the effectiveness of SOMP, we first design ideal frequency-dependent measurement matrices to align the sparse supports of projections at all subcarriers. Then, feasible phase and delay codebooks that can be implemented by hardware including phase shifters and TTD lines are optimized to approximately realize the ideal matrices.
• By exploiting the properties of projections in SOMP in the angular-frequency domain, we further develop several efficient atom selection techniques to reduce the complexity of the proposed DPP algorithm. Specifically, the number of iterations, atoms, and subcarriers involved in atom selection is reduced significantly without sacrificing much performance.
Simulation results illustrate that the proposed approach can be applied alone to achieve satisfactory performance with very low complexity, or serve as a cheap yet effective initializer to improve the convergence performance and speed of the existing high-performance iterative optimization based DPP algorithm.
The rest of this paper is organized as follows. In Section II, the hybrid massive MIMO system, the wideband THz channel model, and the DPP optimization problem are introduced. Then, the extended SSP algorithm is elaborated in Section III, whose complexity is further reduced using the efficient atom selection techniques presented in Section IV. In Section V, the superiority of the proposed approach is validated through numerical results, and the paper is eventually concluded in Section VI.

II. SYSTEM MODEL AND PROBLEM FORMULATION
In this section, we first introduce the hybrid massive MIMO system under the DPP architecture and the practical wideband THz cluster channel model along with beam split. Then, the DPP optimization problem and its conversion to the matrix decomposition form are elaborated.

A. System and Channel Models
We consider the hybrid massive MIMO system illustrated in Fig. 1, where a base station (BS) with an N T -antenna uniform linear array (ULA) and N RF RF chains serves a user with an N Rantenna ULA. We first consider the single-user scenario while the extension to the multi-user scenario will be investigated in our simulation. The number of data streams is N s . Usually, we [24]. To combat the frequency selectivity of channels, OFDM is adopted where the total system bandwidth f s is evenly divided by K subcarriers.  According to measurement results, THz channels usually have clustered characteristics [30].
Considering half-wavelength antenna spacing, we first define the N-dimensional ULA response vector as then the N R × N T channel matrix at the k-th subcarrier from the BS to the user can be expressed as where N c and N p denote the number of clusters and the number of subpaths in a cluster, η is the frequency of the k-th subcarrier with f c and η = fs K denoting the central frequency and the subcarrier frequency spacing, respectively. Besides, α i,j , τ i,j , ψ T i,j,k , ψ R i,j,k denote the complex gain, the delay, the equivalent angle of departure (AoD) at the BS, and the equivalent angle of arrival (AoA) at the user of the j-th subpath in the i-th cluster, respectively. Notice that, the equivalent path angles are dependent on subcarrier index k due to beam split. Specifically, we have ψ i,j,k = f k fc sin(θ i,j ) at both transceiver sides 2 , where sin(θ i,j ) denotes the actual physical path angle. In the i-th cluster, we have denote the mean AoD, mean AoA, and mean delay, respectively, while △θ T j , △θ R j , △τ j follow the Laplacian distribution [30] with zero mean and standard deviation being σ θ T i , σ θ R i , σ τ i , respectively. The standard deviation here is also known as the angular or delay spread. Without loss of generality, we assume α i,j ∈ CN (0, 1), ∀i, j, where τ max denotes the maximum delay. To effectively compensate for beam split, the DPP architecture is adopted. As illustrated in Fig.   1, a small-scale TTD line network is inserted between RF chains and the large-scale phase shifter network. To reduce the hardware cost and energy consumption, each RF chain only connects to where x k denotes the N s ×1 transmit signal at the BS at the k-th subcarrier satisfying E{x k x H k } = 1 Ns I Ns , and ρ denotes the transmit power at the k-th subcarrier 3 . The analog and digital precoders at the k-th subcarrier are denoted by A k ∈ C N T ×N RF and D k ∈ C N RF ×Ns , respectively. The additive white Gaussian noise (AWGN) at the k-th subcarrier is denoted by n k , whose elements follow CN (0, σ 2 n ) with σ 2 n denoting the noise variance. The signal-to-noise ratio (SNR) is defined as ρ/σ 2 n .
The main difference between the DPP architecture and the conventional phase-shifter-based hybrid precoding architecture is that the analog precoder at each subcarrier is frequency-dependent.
Specifically, the analog precoder at the k-th subcarrier can be expressed as where A denotes the common phase matrix introduced by the phase shifter network, and the N T T D × N RF matrix T k denotes the frequency-dependent phases introduced by the TTD line network at the k-th subcarrier. The Kronecker product of T k and 1 M comes from the one-to-M connection mode between TTD lines and antennas at the BS. Denoting the delay matrix of the TTD line network as T , a certain TTD line's phase at a certain subcarrier is proportional to both the delay parameter and the subcarrier frequency and can be expressed as

B. Problem Formulation
We first assume that perfect downlink CSI is available at the BS 4 and the impact of imperfect CSI will be investigated in simulation. Based on the downlink CSI, analog and digital precoders need to be designed to maximize the sum rate of all subcarriers, which can be expressed as In the presence of constraints caused by the hardware used for analog precoding implementation, directly maximizing (6) is difficult due to its complicated form. Therefore, we turn to solve the following matrix decomposition problem to approximate the fully-digital precoders with the hybrid precoders at all subcarriers, which is proved to be approximately equivalent to sum rate maximization while much simpler to solve [10]: In P1, (8) is the unit modulus constraint of phase shifters' phase shifts, t max in (9) denotes the maximal delay of TTD lines, and (10) is to meet the transmit power constraint. To obtain where the N s × 1 power allocation vector, p k , satisfying p k 2 2 = N s , can be obtained by the wellknown water-filling algorithm. For neatness, detailed mathematical procedures of water-filling based power allocation are omitted in this paper.

III. EXTENDED SSP BASED WIDEBAND DPP
In this section, we first briefly introduce the principles of the conventional SSP algorithm designed for narrowband phase-shifter-based hybrid precoding. Then, we extend it properly to solve the considered wideband DPP problem.

A. SSP Based Narrowband Hybrid Precoding
Considering narrowband systems with only phase shifters, P1 is simplified to the following optimization problem with N T T D = 0 and the subscript k discarded: As one of the most classical algorithms to solve P2, the key idea of SSP is that the fully-digital precoder can be approximated by the linear combinations of several array response vectors when the channel is sparse in the angular domain [10], e.g., in the mmWave or THz band. Based on such precoder structures, P2 can be converted to the following CS form: where the N T × G measurement matrix consisting of G atoms (e.g., columns) in the form of ULA response vectors is defined as with the corresponding angles evenly dividing the entire angular space such that Now, let us describe P3 using CS terms. The columns of F are measurements, A is the common measurement matrix, and the columns of the G × N s auxiliary matrix D are sparse vectors corresponding to N s streams that have common sparse supports. The number of non-zero rows in D equals the number of RF chains, which means that the analog precoding vector of each RF chain has to be one of the atoms of the measurement matrix. As a typical MMV-CS problem, P3 can be readily solved by the slightly modified SOMP algorithm demonstrated in Algorithm 1. In each of the N RF iterations, the atom with maximal summed projections of streams on the residual is selected to update the analog precoder while the digital precoder is updated by least square (LS) based on the temporal analog precoder and the fully-digital precoder. At the end, the power of the digital precoder is normalized to meet the transmit power constraint.
As indicated by [24], [25], directly applying SSP to wideband DPP will result in poor performance, which is not surprising because beam split and the delay optimization of TTD lines are not considered in the original designs. Next, we will introduce the proposed extension of SSP from narrowband to the considered wideband DPP scenario.

B. Frequency-Dependent Measurement Matrices
Without beam split and TTD lines, the extension of SSP from narrowband to wideband would be straightforward. Since the equivalent path angles are the same at all subcarriers, the sparse sup- ports of projections would be aligned along the subcarrier dimension. However, with beam split, such common sparsity structure among subcarriers no longer holds if a frequency-independent measurement matrix is still used. As a result, target atoms that have good performance at all subcarriers cannot be accurately captured by simple averaging. Besides, using same atoms as analog precoders of all subcarriers cannot exploit TTD lines' ability of providing frequencydependent phases. To address this issue, we use frequency-dependent measurement matrices 5 , where the measurement matrix used at the k-th subcarrier is defined as with frequency-dependent angles corresponding to atoms being We remark that the factor f k fc in (20) plays the role of predistortion such that the sparse supports of projections at all subcarriers would be aligned after experiencing the distortion of beam split. We would like to mention that we have also tried to use broader beams [19] as atoms to cover wider angular ranges of clusters and optimize the measurement matrices based on channel samples using dictionary learning [34]. Nevertheless, no performance gain is observed since the energy distribution in the angular domain of cluster channels changes from sample to sample, thus cannot be dynamically matched by any fixed measurement matrices. As a result, although these designs have different sample-wise performance, their average performance on a large number of samples is similar.

C. Feasible Phase and Delay Codebooks
So far, feasibility issues are not considered yet. In practice, a group of analog precoders at all subcarriers are feasible only when specific phase and delay parameters satisfying respective constraints can be found. Therefore, we need to further optimize a feasible phase codebook, A, and a feasible delay codebook T , that can be implemented by actual hardware, such that the corresponding feasible measurement matrices are as close to the ideal measurement matrices as possible. Mathematically, the optimization problem is formulated as follows: P4 : min As a special case of P1 when N s = G and the digital precoders are ignored, most existing algorithms for P1 can be used to solve P4 with slight modifications, such as the AM algorithm proposed in [25]. Nevertheless, we propose a simple yet effective approach for codebook design, which can achieve low matrix approximation error with low complexity. To simplify the solution process, the following two Lemmas will be used, which are proved in Appendices A and B. Let us first consider N T T D = N T , in which case the ideal measurement matrices are exactly achievable and the matrix approximation error can be zero. Since phase shifts of phase shifters are frequency-independent, it is straightforward to use the same phase codebook as the narrowband scenario in (17). Denote the initial phase codebook as A ′ , we set A ′ = A. Then, TTD lines are exploited to compensate for frequency-dependent phase differences between A ′ and A ideal k , ∀k. Considering the g-th atom, n T -th antenna, and k-th subcarrier, we should let where a bias of −f c is added to all subcarriers' frequencies according to Lemma 1 and T ′ (n T , g) denotes the unconstrained delay of the n T -th TTD line. Solving (25), we get the initial delay codebook as However, the delay in (26) is negative when φ g < 0. To meet the non-negative constraint of delays, Lemma 2 is used to obtain the final delay codebook as where β = (N T −1) 2fc equals the maximal propagation delay of the electromagnetic wave across the antenna array [24], [26].
Next, we consider the practical case where N T T D ≪ N T and the ideal measurement matrices can only be realized approximately. To minimize the approximation error, the computed delay of the middle antenna of each subarray is shared by all antennas in the subarray such that where n ′ T T D = (2Mn T T D − M)//2 + 1 denotes the index of the n T T D -th subarray's middle antenna in the entire antenna array. When N T T D = N T , M = 1, (28) reduces to (27).
Once the delay codebook is determined, the initial phase codebook needs to be modified accordingly based on Lemma 1 to compensate for the impact of the previously added bias to all subcarriers' frequencies. With slight notation abuse, we still use A to denote the final phase codebook in wideband DPP, which can be expressed as where △ A ′ (n T T D , g) = e j2πfc T (n T T D ,g) , ∀n T T D , g. Eventually, the feasible measurement matrices can be computed by where T k (n T T D , g) = e −j2πf k T (n T T D ,g) is proportional to the true subcarrier frequency. Using A k , ∀k, P1 can be converted to the following MMV-CS problem P7, where the constraint (32) is based on the fact that the sparse supports at all subcarriers are basically aligned even with beam split.

D. DPP Algorithm
Given the above measurement matrices and codebooks, we extend Algorithm 1 to to the digital precoder [19], thus facilitating the calculation of the water-filling based optimal power allocation vector, p eq,k , in line 19.

IV. LOW-COMPLEXITY EXTENDED SSP BASED WIDEBAND DPP
The E-SSP algorithm proposed in the previous section can still have high complexity when the numbers of RF chains, subcarriers, and atoms are large. To enhance the practicality of the proposed approach, in this section, we reduce the complexity of E-SSP dramatically without sacrificing much performance using three efficient atom selection techniques.

A. Non-Iterative Atom Selection
The complexity of SSP is proportional to the number of iterations. In this subsection, we propose to discard the iterative process and select all atoms at once by exploiting the properties of projections in the angular domain. for k = 1 : 1 : K do 5: 10: for k = 1 : 1 : K do 12: end for 15: end for 16 2 V eq,k (:, 1 : N s )diag(p eq,k ); 20: end for Output: A, T , D k , ∀k; In the original SOMP algorithm, the necessity of iteratively subtracting previously selected atoms' contributions comes from the correlation among different atoms. According to (19), the atoms of the ideal measurement matrices 6 are in the form of ULA response vectors. It can be proved that two ULA response vectors corresponding to different angles are asymptotically orthogonal as the length of vectors tends to infinity [5]. Since the number of antennas is relatively large in practical THz massive MIMO systems, the correlation among atoms would be relatively small. Therefore, it is reasonable to directly select the top-N RF atoms based on the initial projections without requiring multiple iterations. Nevertheless, when the angular distance between two atoms is very small, e.g., smaller than the angular Rayleigh resolution of the array [6], their mutual correlation is still non-negligible. To exclude atoms whose large projection comes from leaked power of neighbouring atoms, we further narrow down candidate atoms to only peaks on the initial projection curve 7 , which can be found by various peak-finding algorithms.
To promote intuitive understanding, one example of the initial projection curve and the atoms selected by different methods is given in Fig. 3. As we can see, directly selecting top-N RF atoms results in the selection of multiple nearby atoms with strong mutual correlation, as indicated by the pluses. In contrast, with the additional peak-finding based filtering, the same atoms as the original iterative method are selected, as indicated by the coincident crosses and circles.

B. Hierarchical Atom Selection and Refinement
Although more atoms can improve the performance of E-SSP thanks to higher angular resolution, the complexity grows as well. To enjoy the benefit of high angular resolution with low complexity, in this subsection, we borrow the hierarchical beam sweeping idea from codebook-based precoding approaches [27] to reduce the number of atoms that require projection computation.
Thanks to the continuity of the angles corresponding to atoms, nearby atoms have similar projections so that the projection curve is smooth in the angular domain. Therefore, the relative positions of peaks on the projection curve rarely change with G. For instance, if the 18-th atom is a peak when G = 128, it is very likely that there is also a peak around the 9-th and the 36-th atom when G = 64 and G = 256, respectively. Therefore, in the first stage, coarsegrained measurement matrices with only G c atoms (G c ≪ G), A c k , ∀k, will be used to select N RF coarse-grained atoms, where the columns of each A c k are evenly sampled from the corresponding fine-grained A k with G atoms. Then, 2G a fine-grained atoms will be added around each selected coarse-grained atom, to obtain new measurement matrices, A f k , ∀k. In the second stage, atom selection procedures are executed once again using A f k , ∀k, to determine the final fine-grained atoms. Using this two-stage hierarchical atom selection and refinement method, the number of atoms that require projection computation is at most (G c + 2N RF G a ) while this number is G if the fine-grained A k , ∀k, are directly used. Define the angular resolution reduction ratio as △G G/G c , we will see that the performance degradation is marginal with a large △G and a small G a under typical system settings in simulation, thus reducing the complexity dramatically.

C. Partial Subcarrier Exploitation
Since the number of subcarriers is usually large in practical wideband THz systems due to huge bandwidth, the complexity of E-SSP would be pretty high if projections at all subcarriers need to be computed. In this subsection, we propose to only exploit partial subcarriers' information for atom selection to reduce the complexity similar to [35] which considers wideband channel estimation.
When multiple sparse vectors have the perfect common sparsity structure, averaging the projections of their measurements can reduce the equivalent noise, thus improving performance.
However, this benefit would saturate when the equivalent noise is low enough and no longer a performance bottleneck. For instance, the bottleneck could be the unrecoverable information loss in the imperfect fully-digital precoders in this case. Observing the CS model in (31), we find that the measurements, i.e., the fully-digital precoders at different subcarriers, are computed based on different subchannels. Therefore, the source of measurement noise is actually channel error, which can be caused by either imperfect estimation or feedback in practice. With lower channel error, fewer subcarriers are required to achieve a low enough equivalent noise level. In the extreme case where perfect CSI is available, only a single subcarrier is sufficient to achieve accurate atom selection with noiseless measurements. Nevertheless, the common sparsity structure among subcarriers is still imperfect due to finite angular resolution of atoms and the approximation error between the ideal and feasible measurement matrices. As as result, multiple measurements are still beneficial with perfect CSI while fewer measurements are required compared to the imperfect CSI case, as will be shown in simulation.
Denote the number of exploited subcarriers as K ′ , which can be evenly sampled from all subcarriers, the complexity of atom selection can be reduced by △K times, where △K K/K ′ denotes the subcarrier reduction ratio. After the atoms are selected based on only K ′ subcarriers' fully-digital precoders, the analog precoders at all K subcarriers are determined while the digital precoders at all K subcarriers can still be calculated in the same way. Notice that this technique is tailored for E-SSP and will lead to severe performance degradation if applied to other DPP algorithms, such as the AM algorithm proposed in [25]. The reason is that without utilizing the highly structured codebooks, the optimized phases and delays that have good performance at a few sampled subcarriers are not necessarily good at the remaining subcarriers.

D. Low-Complexity DPP Algorithm
Overall, with the above three efficient atom selection techniques, the low-complexity extended SSP (LCE-SSP) algorithm for wideband DPP is given in

E. Baseline DPP Algorithms and Complexity Analysis
In this subsection, several major baseline DPP algorithms are introduced and the complexity of different DPP algorithms is analyzed and compared, as listed in Table I ψ c = ψ c + ψ c k ′ ; 7: end for 8: g c = P eakF inder(ψ c , 1 : 1 : G c , N RF ); 9: g f = empty vector; 10: for i = 1 : 1 : N RF do 11: 16: 17: 18: 2 V eq,k (:, 1 : N s )diag(p eq,k ); 25: end for Output: A, T , D k , ∀k; precoders (line 18). Notice that the complexity of codebook optimization is not included since it only occurs once offline for each system setting.
Compared to E-SSP, the number of subcarriers involved in the calculation of fully-digital precoders in LCE-SSP reduces from K to K ′ . Besides, the numbers of iterations, subcarriers, and atoms involved for atom selection in Algorithm 3 reduce from N RF , K, G to 1, K ′ , and at most (G c + 2N RF G a ), respectively.
In the AM algorithm proposed in [25], the phases, delays, and digital precoders are updated alternately when the other two are fixed in each iteration. The number of iterations is denoted by N iter and the update of delay parameters is based on the one-dimensional search among S discrete delay grids. Besides, since AM's performance also depends on the initializer, the initialization algorithm's complexity denoted by O Init also needs to be included in its overall complexity.
In the heuristic DPP algorithm proposed in [24], fully-digital precoders are not needed since the analog precoders are computed based on angles of channel paths with very low complexity.
Nevertheless, this algorithm cannot be directly applied to cluster channels where there are far more subpaths than RF chains. Simply using the mean angles of clusters (if assumed available) leads to poor performance. To include heuristic DPP as a meaningful baseline algorithm and enrich the comparison, we propose to find representative angles for it using the LCE-SSP algorithm with slight modifications. Specifically, we set N s to N R and replace the fully-digital precoders with channels in Algorithm 3. Besides, ideal rather than feasible measurement matrices can be used since we are doing pure signal processing. In this way, N RF representative angles can be found such that the constructed pseudo path channels are close to the original cluster channels. We remark that such an extension of heuristic DPP to cluster channels is not trivial since the proposed algorithm is still used. The complexity of finding angles denoted by O F A also needs to be included in the overall complexity of heuristic DPP.
Algorithm Complexity With typical values of system parameters and algorithm hyperparameters, the complexities of E-SSP and AM are comparable while the complexities of LCE-SSP and heuristic DPP are much lower. Intuitively, the average running time of LCE-SSP (several tens of milliseconds) can be dozens or even hundreds of times shorter than that of AM (several seconds) using the same programming language and CPU.
We further remark that in FDD systems, the proposed E-SSP and LCE-SSP also have much lower feedback overhead [10] than the AM algorithm. Thanks to the predetermined measurement matrices and codebooks, only N RF atom indexes are required to configure the hardware to construct the structured analog precoders at the BS, where each index needs log 2 G bits to represent. In contrast, in the AM algorithm, each of the N RF N T T D delay parameters needs log 2 S bits to represent while the N RF N T phase parameters are continuous (or quantized in practice) between 0 and 2π, thus resulting in a much higher feedback overhead.

V. SIMULATION RESULTS
In this section, extensive numerical results 8 are provided to validate the superiority of the proposed approach in terms of performance, complexity, and robustness. Unless specified, the following system and channel parameters will be used as the default setting in simulation: The achievable data rate per subcarrier and the matrix approximation mean-squared error (MSE) defined as  First of all, we can see that LCE-SSP has almost the same performance as E-SSP in both path and cluster channels, demonstrating the effectiveness of the proposed low-complexity techniques. As shown in Fig. 4(a), in path channels, all DPP algorithms achieve close performance to fully-digital precoding in various SNR regimes thanks to the effective compensation for beam split. Further optimization using AM cannot improve the performance of initialization algorithms since they are already near-optimal. In contrast, according to Fig. 4(b), in cluster channels with a 5 • angular spread, the best DPP algorithm only achieves around 90% of the performance upper bound due to more complicated channel distributions. Thanks to our extension, the heuristic DPP algorithm with angles found by the CS algorithm, i.e., LCE-SSP, performs much better than simply using mean angles of clusters. Nevertheless, it is still inferior to all three matrix approximation based algorithms, especially in high SNR regimes since they are not restricted to physical path angles but exploit the fully-digital precoders with richer information to guide the design of hybrid precoders. In AM, since the analog precoders do not have to be in specific forms, the generated beams would have larger degrees of freedom to flexibly cover each channel sample. Or, understood from another point of view, AM can be seen as a dictionary learning algorithm to find the optimal measurement matrices for E-SSP and LCE-SSP, but in a samplewise manner. Although AM has high potentially achievable performance, its actual performance is also dependent on the quality of the initialization point. The role of a good initializer in AM is to construct initial analog precoders that are close to the optimal analog precoders to simplify the iterative optimization process and reduce the probability of falling into bad local optima.
We can see from Fig. 4(b) that LCE-SSP beats the heuristic DPP initializer and the random initializer, therefore can be used as a cheap yet effective initializer for AM in practical cluster THz channels. In Fig. 5, the impact of angular spread is examined at a finer granularity. The performance of E-SSP and LCE-SSP decreases with angular spread since the fully-digital precoders for complicated channels with larger angular spreads are less likely to be well approximated by only a few atoms in the predetermined measurement matrices. Nevertheless, they are still indispensable to guarantee the high performance of AM as good initializers. Besides, when the angular spread is small, e.g., σ θ < 1 • , LCE-SSP alone outperforms AM with the random initializer. The simultaneous good performance and low complexity makes LCE-SSP very appealing to such channel conditions. At last, without exploiting the fully-digital precoders, heuristic DPP is always inferior to LCE-SSP even with CS angles, either applied alone or as AM's initializer. To promote understanding of the importance of beam split compensation, Fig. 6 illustrates the achievable rates with different remaining beam split levels, which are determined by both the number of TTD lines and the fractional bandwidth defined as f s /f c . On the one hand, with larger fractional bandwidth, the initial beam split level is higher, as indicated by different heights of the starting points of the curves. On the other hand, with more TTD lines, beam split can be compensated better, as indicated by the upward trend of the curves. The performance will not improve further when beam split is compensated thoroughly while the minimal number of TTD lines required to achieve this is proportional to the fractional bandwidth [19]. For instance, the performance of LCE-SSP saturates with 32 TTD lines when the fractional bandwidth is 0.1 while only 8 TTD lines are required when the fractional bandwidth is 0.025. In contrast, the curve is still rising at N T T D = 32 when the fractional bandwidth is 0.4.
Next, we show the impact of initializers on AM's convergence performance and speed in Fig.   7. According to Fig. 7(a), both the convergence performance and speed of AM initialized by LCE-SSP are better than that initialized by the other two initializers in two different system scales, indicating its superiority. Same conclusions can be drawn from Fig. 7(b), where the early exit strategy is used to accelerate convergence 9 . Define the exit threshold of matrix approximation MSE in AM as ǫ, we can see that ǫ = 10 −6 is a good choice to decrease the number of iterations  (11) is smaller than that with the heuristic DPP initializer (15) or with the random initializer (18).   respectively. In subfigure (b), the solid curves and the dashed curves denote convergence performance and speed, respectively.
So far we have been assuming perfect CSI is available. However, channel estimation or feedback error is inevitable in practice. So, we investigate the impact of channel error, which is defined as E 1 withĤ k denoting the k-th estimated subchannel, in Fig.   8. From Fig. 8(a), the performance degradation speed of heuristic DPP and LCE-SSP is slower than AM as channel error increases thanks to the restriction on analog precoders' structures.
Specifically, the performance of AM starts to drop quickly when channel error is larger than −15 dB while the performance of the other two algorithms keeps almost unchanged until the channel error is larger than −5 dB. When the channel error is 0 dB, LCE-SSP has similar performance as AM with the heuristic DPP initializer. Such strong robustness to imperfect CSI makes LCE-SSP appealing in challenging scenarios with large noises and limited resources during channel estimation and feedback. Initialized by LCE-SSP, AM still has the best performance thanks to the robust initialization point. As elaborated in Section IV.C, the selection of K ′ depends on the channel error level. From Fig. 8   The extension from the single-user scenario to the multi-user scenario is straightforward in matrix decomposition based DPP algorithms, as long as the fully-digital precoders are available [25]. For the calculation of the fully-digital precoders and the digital precoders at the end, we perform zero-forcing (ZF) with water-filling based power allocation on the original channels and the equivalent channels, respectively. For heuristic DPP, N RF /N U representative angles are contributed by each user's channels to make up all N RF angles, where N U denotes the number of users. We consider each user has a single antenna and one data stream as in [25]. Two RF chain configurations are investigated, where N RF is set to N U to guarantee that each user is served in the first configuration and set to N U N c to thoroughly exploit the spatial multiplexing gain in the second configuration. The performance when N U = 4 is illustrated in Fig. 9. We can see that among all candidate initializers, LCE-SSP always leads to the best performance of AM.
Specifically, it achieves 84% performance of fully-digital precoding in the first configuration and this ratio increases to 95% in the second configuration thanks to sufficient RF chains. Apart from better performance, the LCE-SSP initializer also leads to faster convergence of AM. For instance, when N RF = N U , ǫ = 10 −6 , and SNR = 10 dB, the average number of iterations with the LCE-SSP initializer is about half of that with the random initializer, indicating twice the speedup. In practice, the proposed LCE-SSP can be flexibly applied in different modes according to different requirements and RF chain configurations. When high performance is required and RF chains are limited, LCE-SSP can serve as a cheap yet effective initializer of AM. When low complexity is required and RF chains are sufficient, LCE-SSP alone can already achieve satisfactory performance. We also notice that the advantage of LCE-SSP over heuristic DPP is smaller than the single-user scenario since ZF is not the optimal fully-digital precoder in the multi-user scenario. Nevertheless, the good performance of heuristic DPP is still based on the representative angles found by LCE-SSP.

VI. CONCLUSION
In this paper, we have proposed an extended SSP based DPP algorithm for THz massive MIMO systems. By exploiting the precoder structures, the matrix decomposition problem is converted to an MMV-CS problem, which can be readily solved by the slightly modified SOMP algorithm. To compensate for beam split, ideal frequency-dependent measurement matrices are designed for the alignment of sparse supports at all subcarriers, which can be approximately realized by optimized feasible phase and delay codebooks. To further reduce complexity, several efficient atom selection techniques are developed without sacrificing much performance. According to simulation results, by applied alone or in combination with existing DPP approaches, the proposed approach demonstrates superiority in performance, complexity, and robustness over existing approaches, thus is promising in future 6G THz massive MIMO systems.

APPENDIX
A. Proof of Lemma 1 Consider a TTD line connected to an RF chain has delay t, and the m-th phase shifter connected to it has phase shift e jθm , ∀m, then the superimposed phase of the TTD line and the m-th phase shifter at the k-th subcarrier is e j(θm−2πf k t) , ∀m, k. If we add a bias of −△f to all subcarriers' frequencies, the TTD line's phase at the k-th subcarrier will change from e −j2πf k t to e −j2π(f k −△f )t . To keep the superimposed phases unchanged, we need to modify the phase shifters' phase shifts to e j(θm−2π△f t) , ∀m. The reasoning process holds for any valid △f , thus proving Lemma 1.

B. Proof of Lemma 2
Let us first consider a specific subcarrier and discard the subscript k. Denote the N RF × 1 vector of delay biases added to all TTD lines connected to N RF RF chains as △t, and denote the original analog precoder and the analog precoder with delay biases added as A and A ′ , respectively, we have A ′ = AΣ and A = A ′ Σ −1 , where Σ = diag e −j2πf △t is an invertible diagonal matrix. Based on the properties of Moore-Penrose inverse such that AA † A = A and (AB) † = B † A † , we have Right-multiply both sides of (34) by matrix ΣA ′ † , we have where the property A † AA † = A † , is used to obtain (37). Eventually, based on ΣA ′ † = A † , we prove that AA † = A ′ A ′ † . Since AA † = A ′ A ′ † holds for any valid A, for the k-th subcarrier, i.e., A k D k does not change with △t. Therefore, the sum rate R = K k=1 log 2 I N R + ρ