Low-Complexity Geometry-Based MIMO Channel Simulation

The simulation of electromagnetic wave propagation in time-variant wideband multiple-input multiple-output mobile radio channels using a geometry-based channel model (GCM) is computationally expensive. Due to multipath propagation, a large number of complex exponentials must be evaluated and summed up. We present a low-complexity algorithm for the implementation of a GCM on a hardware channel simulator. Our algorithm takes advantage of the limited numerical precision of the channel simulator by using a truncated subspace representation of the channel transfer function based on multidimensional discrete pro-late spheroidal (DPS) sequences. The DPS subspace representation o ﬀ ers two advantages. Firstly, only a small subspace dimension is required to achieve the numerical accuracy of the hardware channel simulator. Secondly, the computational complexity of the subspace representation is independent of the number of multipath components (MPCs). Moreover, we present an algorithm for the projection of each MPC onto the DPS subspace in O (1) operations. Thus the computational complexity of the DPS subspace algorithm compared to a conventional implementation is reduced by more than one order of magnitude on a hardware channel simulator with 14-bit precision.


INTRODUCTION
In mobile radio channels, electromagnetic waves propagate from the transmitter to the receiver via multiple paths.A geometry-based channel model (GCM) assumes that every multipath component (MPC) can be modeled as a plane wave, mathematically represented by a complex exponential function.The computer simulation of time-variant wideband multiple-input multiple-output (MIMO) channels based on a GCM is computationally expensive, since a large number of complex exponential functions must be evaluated and summed up.
This paper presents a novel low-complexity algorithm for the computation of a GCM on hardware channel simulators.Hardware channel simulators [1][2][3][4][5] allow one to simulate mobile radio channels in real time.They consist of a powerful baseband signal processing unit and radio frequency frontends for input and output.In the baseband processing unit, two basic operations are performed.Firstly, the channel impulse response is calculated according to the GCM.Secondly, the transmit signal is convolved with the channel im-pulse response.The processing power of the baseband unit limits the number of MPCs that can be calculated and hence the model accuracy.We note that the accuracy of the channel simulator is limited by the arithmetic precision of the baseband unit as well as the resolution of the analog/digital converters.On the ARC SmartSim channel simulator [2], for example, the baseband processing hardware uses 16-bit fixedpoint processors and an analog/digital converter with 14-bit precision.This corresponds to a maximum achievable accuracy of E max = 2 −13 .
The new simulation algorithm presented in this paper takes advantage of the limited numerical accuracy of hardware channel simulators by using a truncated basis expansion of the channel transfer function.The basis expansion is based on the fact that wireless fading channels are highly oversampled.Index-limited snapshots of the sampled fading process span a subspace of small dimension.The same subspace is also spanned by index-limited discrete prolate spheroidal (DPS) sequences [6].In this paper, we show that the projection of the channel transfer function onto the DPS subspace can be calculated approximately but very efficiently in O(1) operations from the MPC parameters given by the model.Furthermore, the subspace representation is independent of the number of MPCs.Thus, in the hardware simulation of wireless communication channels, the number of paths can be increased and more realistic models can be computed.By adjusting the dimension of the subspace, the approximation error can be made smaller than the numerical precision given by the hardware, allowing one to trade accuracy for efficiency.Using multidimensional DPS sequences, the DPS subspace representation can also be extended to simulate time-variant wideband MIMO channel models.
One particular application of the new algorithm is the simulation of Rayleigh fading processes using Clarke's [7] channel model.Clarke's model for time-variant frequencyflat single-input single-output (SISO) channels assumes that the angles of arrival (AoAs) of the MPCs are uniformly distributed.Jakes [8] proposed a simplified version of this model by assuming that the number of MPCs is a multiple of four and that the AoAs are spaced equidistantly.Jakes' model reduces the computational complexity of Clarke's model by a factor of four by exploiting the symmetry of the AoA distribution.However, the second-order statistics of Jakes' simplification do not match the ones of Clarke's model [9] and Jakes' model is not wide-sense stationary [10].Attempts to improve the second-order statistics while keeping the reduced complexity of Jakes' model are reported in [6,[9][10][11][12][13][14].However, due to the equidistant spacing of the AoAs, none of these models achieves all the desirable statistical properties of Clarke's reference model [15].Our new approach presented in this paper allows us to reduce the complexity of Clarke's original model by more than an order of magnitude without imposing any restrictions on the AoAs.

Contributions of the paper
(i) We apply the DPS subspace representation to derive a low-complexity algorithm for the computation of the GCM.(ii) We introduce approximate DPS wave functions to calculate the projection onto the subspace in O(1) operations.(iii) We provide a detailed error and complexity analysis that allows us to trade efficiency for accuracy.(iv) We extend the DPS subspace projection to multiple dimensions and describe a novel way to calculate multidimensional DPS sequences using the Kronecker product formalism.
Notation.Let Z, R, and C denote the set of integers, real and complex numbers, respectively.Vectors are denoted by v and matrices by V. Their elements are denoted by v i and V i,l , respectively.Transposition of a vector or a matrix is indicated by • T and conjugate transposition by  [16].
For an N-dimensional, finite index set I ⊂ Z N , the elements of the sequence v m , m ∈ I, may be collected in a vector v.For a parameterizable function f , { f } denotes the family of functions over the whole parameter space.The absolute value, the phase, the real part, and the imaginary part of a complex variable a are denoted by |a|, Φ(a), a, and a, respectively.E {•} denotes the expectation operator.

Organization of the paper
In Section 2, a subspace representation of time-variant frequency-flat SISO channels based on one-dimensional DPS sequences is derived.The main result of the paper, that is, the low-complexity calculation of the basis coefficients of the DPS subspace representation, is given in Section 3. Section 4 extends the DPS subspace representation to higher dimensions, enabling the computer simulation of wideband MIMO channels.A summary and conclusions are given in Section 5. Appendix A proposes a novel way to calculate the multidimensional DPS sequences utilizing the Kronecker product.Appendix B gives a detailed proof of a central theorem.A list of symbols is defined in Appendix C.

Time-variant frequency-flat SISO geometry-based channel model
We start deriving the DPS subspace representation for the generic GCM for time-variant frequency-flat SISO channels depicted in Figure 1.The GCM assumes that the channel transfer function h(t) can be written as a superposition of P MPCs: where each MPC is characterized by its complex weight η p , which embodies the gain and the phase shift, as well as its Florian Kaltenberger et al.Doppler shift ω p .With 1/T S denoting the sampling rate of the system, the sampled channel transfer function can be written as where ν p = ω p T S is the normalized Doppler shift of the pth MPC.We refer to (2) as the sum of complex exponentials (SoCE) algorithm for computing the channel transfer function h m .We assume that the normalized Doppler shifts ν p are bounded by the maximum (one-sided) normalized Doppler bandwidth ν Dmax , which is given by the maximum speed v max of the transmitter, the carrier frequency f C , the speed of light c, and the sampling rate 1/T S , In typical wireless communication systems, the maximum normalized Doppler bandwidth 2ν Dmax is much smaller than the available normalized channel bandwidth (see Figure 2): Thus, the channel transfer function ( 1) is highly oversampled.
Clarke's model [17] is a special case of (2) and assumes that the AoAs ψ p of the impinging MPCs are distributed uniformly on the interval [−π, π) and that E {|η p | 2 } = 1/P.The normalized Doppler shift ν p of the pth MPC is related to the AoA ψ p by ν p = ν Dmax cos(ψ p ). Jakes' model [8] and its variants [9][10][11][12][13][14] assume that the AoAs ψ p are spaced equidistantly with some (random) offset ϑ: If P is a multiple of four, symmetries can be utilized and only P/4 sinusoids have to be evaluated [8].However, the second-order statistics of such models do not match the ones of Clarke's original model [9].In this paper, a truncated subspace representation is used to reduce the complexity of the GCM (2).The subspace representation does not require special assumptions on the AoAs ψ p .It is based on DPS sequences, which are introduced in the following section.

DPS sequences
In this section, one-dimensional DPS sequences are reviewed.They were introduced in 1978 by Slepian [17].Their applications include spectrum estimation [18], approximation, and prediction of band-limited signals [15,17] as well as channel estimation in wireless communication systems [6].DPS sequences can be generalized to multiple dimensions [19].Multidimensional DPS sequences are reviewed in Section 4.2, where they are used for wideband MIMO channel simulation.
Definition 1.The one-dimensional discrete prolate spheroidal (DPS) sequences v (d)  m (W, I) with band-limit W = [−ν Dmax , ν Dmax ] and concentration region I = {M 0 , . . ., M 0 + M − 1} are defined as the real solutions of They are sorted such that their eigenvalues λ d (W, I) are in descending order: To ease notation, we drop the explicit dependence of v (d) m (W, I) on W and I when it is clear from the context.Further, we define the DPS vector v (d) (W, I) ∈ C M as the DPS sequence v (d)  m (W, I) index-limited to I. The DPS vectors v (d) (W, I) are also eigenvectors of the M × M matrix K with elements K m,n = sin(2πν Dmax (m − n))/ π(n − m).The eigenvalues of this matrix decay exponentially and thus render numerical calculation difficult.Fortunately, there exists a tridiagonal matrix commuting with K, which enables fast and numerically stable calculation of DPS sequences [17,20].Figures 3 and 4 illustrate one-dimensional DPS sequences and their eigenvalues, respectively.Some properties of DPS sequences are summarized in the following theorem.

Theorem 1. (1)
The sequences v (d)  m (W, I) are band-limited to W.
(2) The eigenvalue λ d (W, I) of the DPS sequence v (d) m (W, I) denotes the energy concentration of the sequence within I: (5) Every band-limited sequence h m can be decomposed uniquely as h m = h m + g m , where h m is a linear combination of DPS sequences v (d)  m (W, I) for some I and g m = 0 for all m ∈ I.  m , v (1)  m , and v (2)  m for M 0 = 0, M = 256, and Mν Dmax = 2. Proof.See Slepian [17].

DPS subspace representation
The time-variant fading process {h m } given by the model in (2) obtained by index limiting h m to I can be represented as a linear combination of the DPS vectors Properties ( 2) and (3) of Theorem 1 show that the first D = 2ν Dmax M + 1 DPS sequences contain almost all of their energy in the index-set I. Therefore, the vectors {h} span a subspace with essential dimension [6] Due to (4), the time-variant fading process is highly oversampled.Thus the maximum number of subspace dimensions M is reduced by 2ν Dmax 1.In typical wireless communication systems, the essential subspace dimension D is in the order of two to five only.This fact is exploited in the following definition.Definition 2. Let h be a vector obtained by index limiting a band-limited process with band-limit W to the index set I. Further, collect the first D DPS vectors v (d) (W, I) in the matrix The DPS subspace representation of h with dimension D is defined as where α is the projection of the vector h onto the columns of V: For the purpose of channel simulation, it is possible to use D > D DPS vectors in order to increase the numerical accuracy of the subspace representation.The subspace dimension D has to be chosen such that the bias of the subspace representation is small compared to the machine precision of the underlying simulation hardware.This is illustrated in Section 3.2 by numerical examples.
In terms of complexity, the problem of computing the series (2) was reformulated into the problem of computing the basis coefficients α of the subspace representation (13).If they were computed directly using (14), the complexity of the problem would not be reduced.In the following section, we derive a novel low-complexity method to calculate the basis coefficients α approximately.

Approximate calculation of the basis coefficients
In this section, an approximate method to calculate the basis coefficients α in ( 13) with low complexity is presented.Until now we have only considered the time domain of the channel and assumed that the band limiting region W is symmetric around the origin.To make the methods in this section also applicable to the frequency domain and the spatial domains (cf.Section 4), we make the more general assumption that The projection of a single complex exponential vector e p = [e 2π jνpM0 , . . ., e 2π jνp(M0+M−1) ] T onto the basis functions v (d) (W, I) can be written as a function of the Doppler shift ν p , the band-limit region W, and the index set I, Since h can be written as the basis coefficients α (14) can be calculated by where γ p = [γ 0 (ν p ; W, I), . . ., γ D−1 (ν p ; W, I)] T denote the basis coefficients for a single MPC.
To calculate the basis coefficients γ d (ν p ; W, I), we take advantage of the DPS wave functions U d ( f ; W, I).For the special case W 0 = 0 and M 0 = 0 the DPS wave functions are defined in [17].For the more general case, the DPS wave functions are defined as the eigenfunctions of They are normalized such that The DPS wave functions are closely related to the DPS sequences.It can be shown that the amplitude spectrum of a DPS sequence limited to m ∈ I is a scaled version of the associated DPS wave function (cf.[17, equation (26)]) Comparing ( 16) with (21) shows that the basis coefficients can be calculated according to The following definition and theorem show that U d (ν p ; W, I) can be approximately calculated from v (d) m (W, I) by a simple scaling and shifting operation [21].
m (W, I) be the DPS sequences with bandlimit region W = [W 0 − W max , W 0 + W max ] and index set I = {M 0 , . . ., M 0 + M − 1}.Further denote by λ d (W, I) the corresponding eigenvalues.For ν p ∈ W define the index m p by Approximate DPS wave functions are defined as where the sign is taken such that the following normalization holds: Theorem 2. Let ψ d (c, f ) be the prolate spheroidal wave functions [22].Let c > 0 be given and set In other words, both the approximate DPS wave functions as well as the DPS wave functions themselves converge to the prolate spheroidal wave functions.
Proof.For W 0 = 0 and The general case follows by using the two identities Theorem 2 suggests that the approximate DPS wave functions can be used as an approximation to the DPS wave functions.Therefore, the basis coefficients ( 22) can be calculated approximately by The theorem does not indicate the quality of the approximation.It can only be deduced that the approximation improves as the bandwidth W max decreases, while the number of samples M = c/πW max increases.This fact is exploited in the following definition.Definition 4. Let h be a vector obtained by index limiting a band-limited process of the form (2) with band For a positive integer r-the resolution factor-define The approximate DPS subspace representation with dimension D and resolution factor r is given by whose approximate basis coefficients are Note that the DPS sequences are required in a higher resolution only for the calculation of the approximate basis coefficients.The resulting h D,r has the same sample rate for any choice of r.

Bias of the subspace representation
In this subsection, the square bias of the subspace representation bias 2 and the square bias of the approximate subspace representation are analyzed.
For ease of notation, we assume again that W = [−ν Dmax , ν Dmax ], that is, we set W 0 = 0 and W max = ν Dmax .However, the results also hold for the general case (15).If the Doppler shifts ν p , p = 0, . . ., P − 1, are distributed independently and uniformly on W, the DPS subspace representation h coincides with the Karhunen-Loève transform of h [23] and it can be shown that bias 2 If the Doppler shifts ν p , p = 0, . . ., P − 1, are not distributed uniformly, (35) can still be used as an approximation for the square bias [21].
For the square bias of the approximate DPS subspace representation h D,r , no analytical results are available.However, for the minimum achievable square bias, we conjecture that bias 2  min,r = min This conjecture is substantiated by numerical Monte-Carlo simulations using the parameters from Table 1.The Doppler shifts ν p , p = 0, . . ., P − 1, are distributed independently and uniformly on W. The results are illustrated in Figure 5.It can be seen that the square bias of the subspace representation bias 2 h D decays with the subspace dimension.For D ≥ 2Mν Dmax + 1 = 2 this decay is even exponential.These two properties can also be seen directly from (35) and the exponential decay of the eigenvalues λ d (W, I).The square bias bias 2 h D,r of the approximate subspace representation is similar to bias 2 h D up to a certain subspace dimension.Thereafter, the square bias of the approximate subspace representation levels out at bias 2 min,r ≈ (2ν Dmax /r) 2 .Increasing the resolution factor pushes the levels further down.
Let the maximal allowable square error of the simulation be denoted by E 2 max .Then, the approximate subspace representation can be used without loss of accuracy if D and r are chosen such that bias 2 Good approximations for D and r can be found by The first expression can be computed using (35).Using conjecture (36), the latter evaluates to Using a 14-bit fixed-point processor, the maximum achievable accuracy is E 2 max = (2 −13 ) 2 ≈ 1.5 × 10 −8 .For the example of Figure 5, where the maximum Doppler shift ν Dmax = 4.82 × 10 −5 and the number of samples M = 2560, the choice D = 4 and r = 2 makes the simulation as accurate as possible on this hardware.Depending on the application, a lower accuracy might also be sufficient.

Complexity and memory requirements
In this subsection, the computational complexity of the approximate subspace representation (31) is compared to the SoCE algorithm (2).The complexity is expressed in number of complex multiplications (CM) and evaluations of the complex exponential (CE).Additionally, we compare the number of memory access (MA) operations, which gives a better complexity comparison than the actual memory requirements.
We assume that all complex numbers are represented using their real and imaginary part.A CM thus requires four multiplication and two addition operations.As a reference for a CE we use a table look-up implementation with linear interpolation for values between table elements [2].This implementation needs six addition, four multiplication, and two memory access operations.
Let the number of operations that are needed to evaluate h and h be denoted by C h and C h , respectively.Using the SoCE algorithm, for every m ∈ I = {M 0 , . . ., M 0 +M−1} and every p = 0, . . ., P − 1, a CE and a CM have to be evaluated, that is, For the approximate DPS subspace representation with dimension D, first the approximate basis coefficients α have to be evaluated, requiring  operations where the first term accounts for ( 29) and the second term for (32).In total, for the evaluation of the approximate subspace representation (31), operations are required.For large P, the approximate DPS subspace representation reduces the number of arithmetic operations compared to the SoCE algorithm by The memory requirements of the DPS subspace representation are determined by the block length M, the subspace dimension D and the resolution factor r. If the DPS sequences are stored with 16-bit precision, are needed.
In Figure 6, C h and C h are plotted over the number of paths P for the parameters given in Table 1.Multiplications and additions are counted as one operation.Memory access operations are counted separately.The subspace dimension is chosen to be D = 4 according to the observations of the last subsection.The memory requirements for the DPS subspace representation are Mem h = 80 kbyte.
It can be seen that the complexity of the approximate DPS subspace representation in terms of number of arithmetic operations as well as memory access operations increases with slope D, while the complexity of the SoCE algorithm increases with slope M. Since in the given example D M, the approximate DPS subspace representation already enables a complexity reduction by more than one order of magnitude compared to the SoCE algorithm for P = 30 paths.Asymptotically, the number of arithmetic operations can be reduced by a factor of C h /C h → 465.

The wideband MIMO geometry-based channel model
The time-variant GCM described in Section 2.1 can be extended to describe time-variant wideband MIMO channels.
For simplicity we assume uniform linear arrays (ULA) with omnidirectional antennas.Then the channel can be described by the time-variant wideband MIMO channel transfer function h(t, f , x, y), where t denotes time, f denotes frequency, x the position of the transmit antenna on the ULA, y the position of the receive antenna on the ULA [25].
The GCM assumes that h(t, f , x, y) can be written as a superposition of P MPCs, η p e 2π jωpt e −2π jτp f e 2π j/λ sin ϕpx e −2π j/λ sin ψp y , (45) where every MPC is characterized by its complex weight η p , its Doppler shift ω p , its delay τ p , its angle of departure (AoD) ϕ p , and its AoA ψ p (see Figure 7) and λ is the wavelength.More sophisticated models may also include parameters such as elevation angle, antenna patterns, and polarization.
There exist many models for how to obtain the parameters of the MPCs.They can be categorized as deterministic, geometry-based stochastic, and nongeometrical stochastic models [26].The number of MPCs required depends on the scenario modeled, the system bandwidth, and the number of antennas used.In this paper, we choose the number of MPCs such that the channel is Rayleigh fading, except for the lineof-sight component.
For narrowband frequency-flat systems, approximately P 0 = 40 MPCs are needed to achieve a Rayleigh fading statis-tics [13].If the channel bandwidth is increased, the number of resolvable MPCs increases also.The ITU channel models [27], which are used for bandwidths up to 5 MHz in UMTS systems, specify a power delay profile with up to six delay bins.The I-METRA channel models for the IEEE 802.11n wireless LAN standard [28] are valid for up to 40 MHz and specify a power delay profile with up to 18 delay bins.This requires a total number of MPCs of up to P 1 = 18P 0 = 720.Diffuse scattering can also be modeled using a GCM by increasing the number of MPCs.In theory, diffuse scattering results from the superposition of an infinite number of MPCs [29].However, good approximations can be achieved by using a large but finite number of MPCs [30,31].In MIMO channels, the number of MPCs multiplies by N Tx N Rx , since every antenna sees every scatterer from a different AoA and AoD, respectively.For a 4 × 4 system, the total number of MPCs can thus reach up to We now show that the sampled time-variant wideband MIMO channel transfer function is band-limited in time, frequency, and space.Let F S denote the width of a frequency bin and D S the distance between antennas.The sampled channel transfer function can be described as a fourdimensional sequence h m,q,r,s = h(mT S , qF S , rD S , sD S ), where m denotes discrete time, q denotes discrete frequency, s denotes the index of the transmit antenna, and r denotes the index of the receive antenna. 1 Further, let ν p = ω p T S denote the normalized Doppler shift, θ p = τ p F S the normalized delay, ζ p = sin(ϕ p )D S /λ and ξ p = sin(ψ p )D S /λ the normalized angles of departure and arrival, respectively.If all these indices are collected in the vectors m = [m, q, s, r] T , h m can be written as that is, the multidimensional form of (2).The band-limitation of h m in time, frequency, and space is defined by the following physical parameters of the channel.
(1) The maximum normalized Doppler shift of the channel ν Dmax defines the band-limitation in the time domain.It is determined by the maximum speed of the user v max , the carrier frequency f C , the speed of light c, and the sampling rate 1/T S , that is, (2) The maximum normalized delay of the scenario θ max defines the band-limitation in the frequency domain.
It is determined by the maximum delay τ max and the sample rate 1/F S in frequency (3) The minimum and maximum normalized AoA, ξ min and ξ max define the band-limitation in the spatial domain at the receiver.They are given by the minimum and maximum AoA, ψ min and ψ max , the spatial sampling distance D S and the wavelength λ: The band-limitation at the transmitter is given similarly by the normalized minimum and maximum normalized AoD, ζ min and ζ max .
In summary it can be seen that h m is band-limited to Thus the discrete time Fourier transform (DTFT) vanishes outside the region W, that is,

Multidimensional DPS sequences
The fact that h m is band-limited allows one to extend the concepts of the DPS subspace representation also to time-variant wideband MIMO channels.Therefore, a generalization of the one-dimensional DPS sequences to multiple dimensions is required. where They are sorted such that their eigenvalues λ d (W, I) are in descending order To ease notation, we drop the explicit dependence of v (d) m (W, I) on W and I when it is clear from the context.Further, we define the multidimensional DPS vector v (d) (W, I) ∈ C L as the multidimensional DPS sequence v (d) m (W, I) index-limited to I. In particular, if every element m ∈ I is indexed lexicographically, such that I = {m l , l = 0, 1, . . ., L − 1}, then All the properties of Theorem 1 also apply to multidimensional DPS sequences [19].The only difference is that m has to be replaced with m and Z with Z N .
Example 1.In the two-dimensional case N = 2 with bandlimiting region W and index set I given by Equation ( 54) reduces to Note that due to the nonsymmetric band-limiting region W, the solutions of (59) can take complex values.Examples of two-dimensional DPS sequences and their eigenvalues are given in Figures 8 and 9, respectively.They have been calculated using the methods described in Appendix A.

Multidimensional DPS subspace representation
We assume that for hardware implementation, h m is calculated blockwise for M samples in time, Q bins in frequency, N Tx transmit antennas, and N Rx receive antennas.Accordingly, the index set is defined by The DPS subspace representation can easily be extended to multiple dimensions.Let h be the vector obtained by index limiting the sequence h m (47) to the index set I (60) and sorting the elements lexicographically.In analogy to the one-dimensional case, the subspace spanned by {h} is also spanned by the multidimensional DPS vectors v (d) (W, I) defined in Section 4.2.Due to the common notation of oneand multidimensional sequences and vectors, the multidimensional DPS subspace representation of h can be defined similarly to Definition 2.  Definition 6.Let h be a vector obtained by index limiting a multidimensional band-limited process of the form (47) with band-limit W to the index set I. Let v (d) (W, I) be the multidimensional DPS vectors for the multidimensional band-limit region W and the multidimensional index set I. Further, collect the first D DPS vectors v (d) (W, I) in the matrix

EURASIP Journal on Advances in Signal Processing
The multidimensional DPS subspace representation of h with subspace dimension D is defined as where α is the projection of the vector h onto the columns of V: The subspace dimension D has to be chosen such that the bias of the subspace representation is small compared to the machine precision of the underlying simulation hardware.The following theorem shows how the multidimensional projection (63) can be reduced to a series of onedimensional projections.

Theorem 3. Let h D be the N-dimensional DPS subspace representation of h with subspace dimension D, band-limiting region W, and index set I. If W and I can be written as Cartesian products
where W i = [W 0,i − W max,i , W 0,i + W max,i ], and I i = {M 0,i , . . ., M 0,i + M i − 1}, then for every d = 0, . . ., D − 1, there exist d 0 , . . ., d N−1 such that the N-dimensional DPS basis vectors v (d) (W, I) can be written as Further, the basis coefficients of the approximate DPS subspace representation are given by where γ (i) p,d = γ di ( f p,i , W i , I i ) are the one-dimensional approximate basis coefficients defined in (29).Additionally, resolution factors r i can be used to improve the approximation.

Proof. See Appendix B
The band-limiting region W (51) and the index set I (60) of the channel model (47) fulfill the prerequisites of Theorem 3 with Thus, Theorem 3 allows us to use the methods of Section 3.1 to calculate the basis coefficients of the multidimensional DPS subspace representation approximately with low complexity.The resolution factors r i , i = 0, . . ., N − 1, have to be chosen such that the bias of the subspace representation is small compared to the machine precision E max of the underlying simulation hardware.A necessary but not sufficient condition for this is to use the methods of Section 3.2 for each dimension independently, that is, to choose r i = 2W max,i /E max .However, it has to be verified numerically that the multidimensional DPS subspace representation achieves the required numerical accuracy.

Complexity and memory requirements
In this subsection, we evaluate the complexity and memory requirements of the N-dimensional SoCE algorithm and the N-dimensional approximate DPS subspace representation, given by Theorem 3.These results are a generalization of the results of Section 3.3.We assume that the one-dimensional DPS sequences v (di) (W i , I i ), i = 0, . . ., N − 1, have been precalculated.Further, we assume that Let the number of operations that are needed to evaluate h (47) and h D (67) be denoted by C h and C h D , respectively.For the SoCE algorithm, For the approximate DPS subspace representation with dimension D, firstly the N-dimensional DPS basis vectors need to be calculated from the one-dimensional DPS vectors (cf.(66)), requiring Secondly, the approximate basis coefficients α have to be evaluated according to (68), requiring In total, for the evaluation of the approximate subspace representation (67), operations are required.Asymptotically for P → ∞, the N-dimensional DPS subspace representation reduces the number of arithmetic operations compared to the SoCE algorithm by the factor The memory requirements of the DPS subspace representation are determined by the size of the index set I, the number of DPS vectors D i , and the resolution factors r i .If the DPS sequences are stored with 16-bit precision, are needed.

Numerical examples
Section 3 demonstrated that an application of the approximate DPS subspace representation to the time-domain of wireless channels may save more than an order of magnitude in complexity.In this subsection, the multidimensional approximate DPS subspace representation is applied to an example of a time-variant frequency-selective channel as well as an example of a time-variant frequency-selective MIMO channel.A comparison of the arithmetic complexity is given.We assume a 14-bit fixed-point hardware architecture, that is, a maximum allowable square error of  1.We assume a typical urban environment with a maximum delay spread of τ max = 3.7 milliseconds given by the ITU Pedestrian B channel model [27].By omitting the spatial domains x and y in (47), we obtain a time-variant frequency-selective GCM where m = [m, q] T and f p = [ν p , θ p ] T .Since ( 76) is bandlimited to and we wish to calculate (76) in the index set we can apply a two-dimensional DPS subspace representation (Definition 6) to (76).Further, we can use Theorem 3 to calculate the basis coefficients α of the subspace representation.
For a given maximum allowable square bias E 2 max = (2 −13 ) 2 , the estimated values of the resolution factors in the time and frequency domain are r 0 = 2ν Dmax /E max ≈ 2 and r 1 = θ max /E max ≈ 512 (rounded to the next power of two).The square bias bias 2 of the two-dimensional exact and the approximate DPS subspace representation is plotted in Figure 10 against the subspace dimension D. It can be seen that bias 2 h D ≈ E 2 max at a subspace dimension of approximately D = 80.The maximum number of one-dimensional DPS vectors is D 0 = 4 and D 1 = 23.

Time, frequency, and spatial domain
Table 3 contains the simulation parameters of the numerical experiments in the spatial domain.The remaining parameters are chosen according to Tables 1 and 2. We assume uniform linear arrays at the transmitter and the receiver with   spacing D S = λ/2 and N Tx = N Rx = 8 antennas each.Further we assume that there is only one cluster of scatterers in the scenario which is not in the vicinity of the transmitter or receiver (see Figure 11) and we assume no line-of-sight component.The AoD and AoA are assumed to be limited by [ϕ min , ϕ max ] = [ψ min , ψ max ] = [−5 • , 5 • ], which has been observed in measurements [32].
A four-dimensional DPS subspace representation is applied to the channel transfer function (47) with W and I defined in (51) and (60).Following the same procedure as in the previous subsection, for a numerical accuracy of 14 bits the estimated values of the resolution factors and the number of one-dimensional DPS vectors in the spatial domains are r 2 = (ζ max − ζ min )/E max ≈ 512, r 3 = (ξ max − ξ min )/E max ≈ 512 (rounded to the next power of 2), and D 2 = D 3 = 5.

Hybrid DPS subspace representation
Last but not least, we propose a hybrid DPS subspace representation that applies a DPS subspace representation in time where the band-limit region W and the index set I are the same as in the two-dimensional case (cf.( 77) and ( 78)).Then, the two-dimensional DPS subspace representation can be applied to each h s,r m , s = 0, . . ., N Tx − 1, r = 0, . . ., N Rx − 1, independently.

Results and discussion
A complexity comparison of the SoCE algorithm and the approximate DPS subspace representation for one, two, and four dimensions is given in Figure 12.It was evaluated using (70) and (73).Also shown is the complexity of the four-dimensional hybrid DPS subspace representation.It can be seen that for time-variant frequency-flat SISO channels, the one-dimensional DPS subspace representation requires fewer arithmetic operations for P > 2 MPCs.The more MPCs are used in the GCM, the more complexity is saved.Asymptotically, the number of arithmetic operations is reduced by C h /C h → 465.
For time-variant frequency-selective SISO channels, the two-dimensional DPS subspace representation requires fewer arithmetic operations for P > 30 MPCs.However, as noted in Section 4.1, channel models for systems with the given parameters require P = 400 paths or more.For such a scenario, the DPS subspace representation saves two orders of magnitude in complexity.Asymptotically, the number of arithmetic operations is reduced by a factor of C h /C h → 6.8 × 10 3 (cf.( 74)).The memory requirements are Mem h = 5.83 Mbyte (cf.(75)).
For time-variant frequency-selective MIMO channels, the four-dimensional DPS subspace representation requires fewer arithmetic operations for P > 2 × 10 3 MPCs.Since MIMO channels require the simulation of up to 10 4 MPCs  (cf.Section 4.1), complexity savings are still possible.The asymptotic complexity savings are C h /C h → 1.9 × 10 4 .However, in the region P < 2 × 10 3 MPCs, the four-dimensional DPS subspace representation requires more complex operations than the corresponding SoCE algorithm.Thus, even though we choose a "best case" scenario with only one cluster, a small angular spread and a low numerical accuracy, there is hardly any additional complexity reduction if the DPS subspace representation is applied in the spatial domain.
The hybrid DPS subspace representation on the other hand exploits the savings of the DPS subspace representation in the time and frequency domain only.From Figure 12 it can be seen that it has fewer arithmetic operations than the four-dimensional DPS subspace representation and the fourdimensional SoCE algorithm for 60 < P < 2 × 10 3 MPCs.Thus the hybrid method is preferable for channel simulations in this region.Further, this method also allows for an efficient parallelization on hardware channel simulators [33].

CONCLUSIONS
We have presented a low-complexity algorithm for the computer simulation of geometry-based MIMO channel models.The algorithm exploits the low-dimensional subspace spanned by multidimensional DPS sequences.By adjusting the dimension of the subspace, it is possible to trade computational complexity for accuracy.Thus the algorithm is ideally suited for fixed-point hardware architectures with limited precision.
We demonstrated that the complexity reduction depends mainly on the normalized bandwidth of the underlying fading process in time, frequency, and space.If the bandwidth is very small compared to the sampling rate, the essential subspace dimension of the process is small and the complexity can be reduced substantially.In the time domain, the maximum Doppler bandwidth of the fading process is much smaller than the system sampling rate.Compared with the SoCE algorithm, our new algorithm reduces the complexity by more than one order of magnitude on 14-bit hardware.
The bandwidth of a frequency-selective fading process is given by the maximum delay in the channel, which is a factor of five to ten smaller than the sampling rate in frequency.Therefore, the DPS subspace representation also reduces the computational complexity when applied in the frequency domain.To achieve a satisfactory numerical accuracy, the resolution factor in the approximation of the basis coefficients needs to be large, resulting in high memory requirements.On the other hand, it was shown that the number of memory access operations is small.Since this figure has more influence on the run-time of the algorithm, the approximate DPS subspace representation is preferable over the SoCE algorithm for a frequency-selective fading-process.
The bandwidth of the fading process in the spatial domain is determined by the angular spread of the channel, which is almost as large as the spatial sampling rate for most scenarios in wireless communications.Therefore, applying the DPS subspace representation in the spatial domain does not achieve any additional complexity reduction for the scenarios of interest.As a consequence, for the purpose of wideband MIMO channel simulation, we propose to use a hybrid method which computes the complex exponentials in the spatial domain directly and applies the subspace representation to the time and frequency domain only.This method also allows for an efficient parallelization on hardware channel simulators.

APPENDICES A. CALCULATION OF MULTIDIMENSIONAL DPS SEQUENCES
In the one-dimensional case (N = 1), where W = [W 0 − W max , W 0 + W max ] and I = {M 0 , . . ., M 0 + M − 1}, the DPS sequences can be calculated efficiently [17,20].The efficient and numerically stable calculation of multidimensional DPS sequences with arbitrary W and I is not trivial and has not been treated satisfactorily in the literature.In this section a new way of calculating multidimensional DPS sequences is derived if their passband region can be written as a Cartesian product of one-dimensional intervals.Indexing every element m ∈ I lexicographically, such that I = {m l , l = 0, 1, . . ., L − 1}, we define the matrix K (W) by where the kernel K (W) is given by (55).Let v (d) (W, I) and λ d (W, I), d = 0, . . ., L−1, denote the eigenvectors and eigenvalues of K (W) : where It can be shown that the eigenvectors v (d) (W, I) and the eigenvalues λ d (W, I) are exactly the multidimensional DPS vectors defined in (57) and their corresponding eigenvalues.
If the DPS sequences are required for m / ∈ I, they can be extended using (54).
The multidimensional DPS vectors can theoretically be calculated for an arbitrary passband region W directly from the eigenproblem (A.2).However, since the matrix K (W)  has an exponentially decaying eigenvalue distribution, this method is numerically unstable.
If W can be written as a Cartesian product of onedimensional intervals (i.e., W is a hyper-cube), where , and the index-set I is written as where I i = {M 0,i , . . ., M 0,i + M i − 1}, the defining kernel K (W)  for the multidimensional DPS vectors evaluates to where u = [u 0 , . . ., u N−1 ] T ∈ I.This means that the kernel K (W) is separable and thus the matrix K (W) can be written as a Kronecker product where K (Wi) , i = 0, . . ., N − 1, are the kernel matrices corresponding to the one-dimensional DPS vectors.Now let λ di (W i , I i ) and v (di) (W i , I i ), d i = 0, . . ., M i − 1, denote the eigenvalues and the eigenvectors of K (Wi) , i = 0, . . ., N − 1, respectively.Then the eigenvalues of K (W) are given by [34,Chapter 9]

B. PROOF OF THEOREM 3
For I given by (65), h can be written as where e (i) p = [e 2π j fp,iM0,i , . . ., e 2π j fp,i(M0,i+Mi−1) ] T .Further, since W is given by (64), the results of Appendix A can be used and V can be written as where every M i × D i matrix V i contains the one-dimensional DPS vectors v d (W i , I i ) in its columns.Using the identity

Figure 2 :
Figure 2: Doppler spectrum H(ν) of the sampled time-variant channel transfer function h m .The maximum normalized Doppler bandwidth 2ν Dmax is much smaller than the available normalized channel bandwidth.

) ( 3 )
The eigenvalues λ d (W, I) satisfy 1 < λ i (W, I) < 0. They are clustered around 1 for d ≤ D − 1, and decay exponentially for d ≥ D , where D = |W ||I| + 1. (4) The DPS sequences v (d) m (W, I) are orthogonal on the index set I and on Z.

Figure 6 :
Figure6: Complexity in terms of number of arithmetic operations (left abscissa) and memory access operations (right abscissa) versus the number of MPCs P. We show results for the sum of complex exponentials algorithm (denoted by "SoCE") and the approximate subspace representation (denoted by "DPSS") using M = 2560, ν Dmax = 4.82 × 10 −5 , and D = 4.

ψ 1 ψ 2 Figure 7 :
Figure 7: Multipath propagation model for a time-variant wideband MIMO radio channel.The signals sent from the transmitter, moving at speed v, arrive at the receiver.Each path p has complex weight η p , time delay τ p , Doppler shift ω p , angle of departure ϕ p , and angle of arrival ψ p .

Definition 5 .
Let I ⊂ Z N be an N-dimensional finite index set with L = |I| elements, and W ⊂ (−1/2, 1/2) N an Ndimensional band-limiting region.Multidimensional discrete prolate spheroidal (DPS) sequences v(d)  m (W, I) are defined as the solutions of the eigenvalue problem

Figure 10 :
Figure 10: bias 2 h D for the subspace representation in the time and frequency domain with ν Dmax = 4.82 × 10 −5 , M = 2560, θ max = 0.056, and Q = 256.The resolution factors are fixed to r 0 = 2 and r 1 = 512.The thin horizontal line denotes the numerical accuracy of a fixed-point 14-bit processor.

Figure 11 :
Figure 11: Scenario of a mobile radio channel with one cluster of scatterers.The AoD and the AoA are limited within the intervals Φ = [ϕ min , ϕ max ] and Ψ = [ψ min , ψ max ], respectively.

Figure 12 :
Figure12: Complexity in terms of number of arithmetic operations versus the number of MPCs P. We show results for the SoCE algorithm (denoted by "SoCE") and the approximate DPS subspace representation (denoted by "DPSS") for one, two, and four dimensions.Also shown is the complexity of the four-dimensional hybrid DPS subspace representation (denoted by "Hybrid").

3 ) 1 e
the basis coefficients α can be calculated byα = V H h = (0) p ⊗ • • • ⊗ e (N−1) p ⊗ • • • ⊗ V H N−1 e (N−1) Figure 1: GCM for a time-variant frequency-flat SISO channel.Signals sent from the transmitter, moving at speed v, arrive at the receiver via different paths.Each MPC p has complex weight η p and Doppler shift ω p • H .The Euclidean ( 2 ) norm of the vector a is denoted by a .The Kronecker product and the Khatri-Rao product (columnwise Kronecker product) are denoted by ⊗ and , respectively.The inner product of two vectors of length N is defined asx, y = N−1 i=0 x i y * i , where • * denotes complex conjugation.If X is a discrete index set, |X| denotes the number of el- is band-limited to the region W = [−ν Dmax , ν Dmax ].Let I = {M 0 , . . ., M 0 + M − 1} denote a finite index set on which we want to calculate h m .Due to property (5) of Theorem 1, h m can be decomposed into h m = h m +g m , where h m is a linear combination of the DPS sequences v (d) m (W, I) and h m = h m for all m ∈ I. Therefore, the vectors

Table 1 :
[24]lation parameters for the numerical experiments in the time domain.The carrier frequency and the sample rate resemble those of a UMTS system[24].The block length is chosen to be as long as a UMTS frame.

Table 2 :
Simulation parameters for the numerical experiments in the frequency domain.Table 2 contains the simulation parameters of the numerical experiments in the frequency domain.The parameters in the time domain are chosen according to Table

Table 3 :
Simulation parameters for the numerical experiments in the spatial domains.
, ν p : Doppler shift and normalized Doppler shift of the pth MPC ω Dmax , ν Dmax : Maximum Doppler shift, maximum normalized Doppler shift τ p , θ p : Delay and normalized delay of the pth MPC τ max , θ max : (ν), U d (ν): DPS wave function and approximate DPS wave function α d , α d : dth basis coefficient and approximate basis coefficient of DPS subspace representation of h γ p,d , γ p,d : dth basis coefficient and approximate basis coefficient of DPS subspace representation of the pth MPC r i , D i :Resolution factor and maximum number of one-dimensional DPS vectors in time (i = 0), frequency (i = 1), space at the transmitter (i = 2), and space at the receiver (i = 3)