Large-Scale MIMO Detection for 3GPP LTE: Algorithms and FPGA Implementations

I Introduction

Multiple-input multiple-output (MIMO) in combination with spatial multiplexing builds the foundation of most modern wireless communication standards, such as 3GPP LTE or IEEE 802.11n. MIMO technology offers significantly higher data rates over single-antenna systems by transmitting multiple data streams concurrently and in the same frequency band. Conventional MIMO wireless systems, however, already start to approach their throughput limits. Consequently, the deployment of novel transceiver technologies is of paramount importance in order to meet the ever-growing demand for higher data rates, better link reliability, and improved coverage, without further increasing the communication bandwidth .

Large-scale (or massive) MIMO is an emerging technology, which postulates the use of antenna arrays having orders of magnitude more elements at the base station (BS) compared to conventional (small-scale) MIMO systems, while serving tens of users simultaneously and in the same frequency band . This technology promises significant improvements in terms of spectral efficiency, link reliability, and coverage compared to conventional (small-scale) systems .

Unfortunately, the promised benefits of large-scale MIMO come at the cost of significantly increased computational complexity in the BS, as opposed to small-scale MIMO systems, which commonly deploy $2$-to-$4$ antennas at both ends of the wireless link. In particular, data detection in the large-scale MIMO uplink is expected to be among the most critical tasks in terms of complexity and power consumption, as the presence of hundreds of antennas at the BS and a large number of users will increase the computational complexity by orders of magnitude. In addition, current cellular systems, such as 3GPP-LTE or LTE-Advanced (LTE-A) , rely on single-carrier frequency division multiple access (SC-FDMA), which further increases the dimensionality (and hence the complexity) of the underlying detection problem. As a consequence, optimal data detection methods, such as maximum-likelihood (ML) detection or soft-output sphere decoding (SD) , whose (average) computational complexity scales exponentially in the number of transmitted data streams , would simply result in prohibitive complexity. Hence, one has to resort to low-complexity (but sub-optimal) linear detection schemes or stochastic detection algorithms that deliver acceptable error-rate performance and scale favorably to the high-dimensional detection problems faced in SC-FDMA-based large-scale MIMO systems.

I-B Contributions

This paper addresses the complexity issue of data detection in SC-FDMA-based large-scale MIMO systems in the uplink, i.e., where multiple users communicate with the BS. We focus on linear soft-output detection in combination with a new approximate matrix inversion method relying on a Neumann series expansion, which significantly reduces the computational complexity compared to that of an exact matrix inversion method. We analyze the implementation trade-offs associated with approximate and exact linear data detection in the large-scale MIMO uplink, and we show analytically that the approximation error caused by the proposed approximate inversion method depends on the ratio between BS antennas and users. We show that the proposed approximation performs well for medium to large ratios between BS antennas and users, while exact linear detection is advantageous for small antenna ratios. We present reference FPGA designs for both, the approximate and exact matrix inversion, and for various antenna configurations, which enables us to characterize the associated hardware complexity vs. error-rate performance trade-offs. The resulting FPGA designs are—to the best of our knowledge—the first data detection engines for massive MIMO systems reported in the open literature that achieve a peak uplink throughput exceeding the $300$ Mb/s specified in 3GPP LTE-Advanced operating at 20 MHz bandwidth .

I-C Notation

I-D Paper Outline

The remainder of the paper is organized as follows. Section II introduces the uplink system model and outlines the basics of linear detection for SC-FDMA-based systems. The approximate matrix inversion approach, a corresponding error analysis, and an error-rate performance/complexity comparison are shown in Section III. Section IV details our VLSI architecture. Section V provides reference FPGA implementation results and a trade-off analysis. We conclude in Section VI. All proofs are relegated to the Appendices.

II Large-Scale MIMO in LTE Uplink

We next introduce the LTE uplink model and present a new and efficient method for linear soft-output minimum mean-square error (MMSE) detection in SC-FDMA-based systems.

We consider the large-scale multi-user (MU) MIMO uplink with ${B}$ antennas at the base-station (BS) communicating with ${U}\leq{B}$ single-antenna users.More generally, instead of having ${U}$ single-antenna users, the proposed system may equivalently support ${U}$ spatial streams, which can, for example, be shared among a smaller number of user terminals that are equipped with more than one antenna. To reduce the peak-to-average power ratio of the user equipment, LTE uplink employs SC-FDMA (short for single-carrier frequency division multiple access) . The ${U}$ users first encode their own transmit bits using channel encoders and then, map the coded bit stream to time-domain constellation points in the finite alphabet $\mathcal{O}$ with cardinality $M=\left|\mathcal{O}\right|$ and average transmit power $E_{s}$ per symbol. An $L$-point discrete Fourier transform (DFT) blockIn practice, the DFT and inverse DFT are carried out by fast (inverse) Fourier transform (I/FFT) units. is used to perform modulation of these time-domain symbols onto orthogonal frequency bands. The $L$ time-domain constellation points for the $i^{\text{th}}$ user are subsumed in the vector \mathbf{x}^{(i)}=\big{[}x^{(i)}_{1},\ldots,x^{(i)}_{L}\big{]}^{T}. The output of the DFT block, namely the frequency-domain symbol, is defined as \mathbf{s}^{(i)}=\big{[}\mathbf{s}^{(i)}_{1},\ldots,\mathbf{s}^{(i)}_{L}\big{]}^{T}=\mathbf{F}_{L}\mathbf{x}^{(i)}. Subsequent processing performed for each user corresponds to that of conventional orthogonal frequency-division multiplexing (OFDM) transmission . Specifically, for each user, the frequency-domain symbols are first mapped onto data-carrying subcarriers and then, transformed back to the time domain with an inverse DFT (IDFT). After prepending the cyclic prefix to the time-domain symbols, all ${U}$ users transmit their time-domain signals simultaneously over the wireless channel.

At the BS, each receive antenna obtains a mixture of the time-domain signals from all users. For data detection, the time-domain signals received at each antenna are first transformed back into the frequency domain using a DFT. The data-carrying symbols are then extracted from the DFT’s output. Assuming a sufficiently long cyclic-prefix (i.e., longer than the delay spread of the channel’s impulse response), the received frequency-domain symbols can be modeled using the standard input-output relation $\mathbf{y}=\mathbf{H}\mathbf{s}+\mathbf{n}$, with the following definitions:

Here, the vector \mathbf{y}^{(i)}=\big{[}y^{(i)}_{1},\ldots,y^{(i)}_{L}\big{]}^{T} contains the received symbols on the $i^{\text{th}}$ antenna in the frequency domain, where $y^{(i)}_{w}$ is the symbol received on the $w^{\text{th}}$ subcarrier of the $i^{\text{th}}$ antenna. The $L\times L$ diagonal matrix \mathbf{H}^{(i,j)}=\text{diag}\big{(}h^{(i,j)}_{1},\ldots,h^{(i,j)}_{L}\big{)} contains the channel’s frequency response of length $L$ between the $i^{\text{th}}$ receive antenna and $j^{\text{th}}$ transmit antenna on its main diagonal, and \mathbf{n}^{(i)}=\big{[}n^{(i)}_{1},\ldots,n^{(i)}_{L}\big{]}^{T} models thermal noise at the $i^{\text{th}}$ receive antenna in the frequency domain. The entries of the vector $\mathbf{n}^{(i)}$ are assumed to be i.i.d. zero-mean Gaussian with variance $N_{0}$ per complex entry.

II-B Linear MMSE Detection

The task of a data detector for MIMO systems is to compute soft-estimates in the form of log-likelihood ratio (LLR) values for each coded bit, given the channel matrixIn practice, channel-state information is acquired using pilot sequences specified by the standard . For the sake of simplicity, we assume perfect channel state information (CSI) throughout the paper. An investigation of the impact of imperfect CSI on the error-rate performance is left for future work. $\mathbf{H}$ and receive vector $\mathbf{y}$. In order to arrive at low computational complexity for data detection in SC-FDMA-based large-scale MIMO systems, we focus exclusively on linear soft-output detection . Linear detection for SC-FDMA mainly consists of the following two steps: (i) channel equalization to generate estimates of the frequency domain symbols, and (ii) soft-output computation to generate LLRs from the equalized frequency domain symbols. Both of these steps are detailed next.

The most common approach to linear MIMO detection is the minimum-mean square error (MMSE) equalizer, which computes equalized frequency-domain symbols as $\hat{\mathbf{s}}=\mathbf{W}\mathbf{y}$ with the MMSE equalization matrix defined as follows :

Since the effective channel matrix $\mathbf{H}$ is built from diagonal $L\times L$ submatrices, we can apply MMSE equalization on a per-subcarrier basis. Specifically, the received frequency symbols on the $w^{\text{th}}$ subcarrier in the frequency domain can be modeled as $\mathbf{y}_{w}=\mathbf{H}_{w}\mathbf{s}_{w}+\mathbf{n}_{w}$, where

Here, $y^{(i)}_{w}$ is the frequency symbol received on the $w^{\text{th}}$ subcarrier for the $i^{\text{th}}$ antenna, and $h^{(i,j)}_{w}$ is the frequency gain (or attenuation) on the $w^{\text{th}}$ subcarrier between the $i^{\text{th}}$ receive antenna and $j^{\text{th}}$ transmit antenna. The scalar $s^{(j)}_{w}$ denotes the symbol transmitted by the $j^{\text{th}}$ user on the $w^{\text{th}}$ subcarrier; the scalar $n^{(i)}_{w}$ models thermal noise at the $i^{\text{th}}$ receive antenna on the $w^{\text{th}}$ subcarrier. With this reformulation, the equalized symbols on the $w^{\text{th}}$ subcarrier are given by $\hat{\mathbf{s}}_{w}=\mathbf{W}_{w}\mathbf{y}_{w}$, with the per-subcarrier MMSE equalization matrix defined as

(ii) LLR computation

To obtain symbol estimates in the time domain, the MMSE detector performs an IDFT on the equalized frequency domain symbols for each user. The time-domain symbol estimates for the $i^{\text{th}}$ user are given by $\hat{\mathbf{x}}^{(i)}=\mathbf{F}_{L}^{H}\hat{\mathbf{s}}^{(i)}$, where $\mathbf{F}_{L}^{H}$ is the IDFT matrix and \hat{\mathbf{x}}^{(i)}=\big{[}\hat{x}^{(i)}_{1},\ldots,\hat{x}^{(i)}_{L}\big{]}^{T} contains the time-domain symbol estimates of the symbols transmitted by the $i^{\text{th}}$ user. To extract LLRs from the time-domain symbol estimates, we approximate each estimate as an independent Gaussian random variable. In particular, the estimated $t^{\text{th}}$ symbol transmitted from the $i^{\text{th}}$ user is modeled as $\hat{x}^{(i)}_{t}=\mu^{(i)}x^{(i)}_{t}+{e^{(i)}_{t}}$, where $\mu^{(i)}$ is the effective channel gain and $e^{(i)}_{t}$ is the post-equalization noise-plus-interference (NPI) variance. Let $\nu_{i}^{2}$ be the variance of $e^{(i)}_{t}$ and $b$ be the bit index of the LLR associated with the $t^{\text{th}}$ symbol transmitted from the $i^{\text{th}}$ user. With this model, the max-log LLRs can be computed as

where $\rho^{2}_{i}={\left(\mu^{(i)}\right)^{2}}/{\nu_{i}^{2}}$ is the post-equalization signal-to-noise-plus-interference ratio (SINR), and $\mathcal{O}_{b}^{0}$ and $\mathcal{O}_{b}^{1}$ correspond to the sets of constellation symbols for which the $b^{\text{th}}$ bit equals to and $1$, respectively.

In order to obtain an explicit formulation of the effective channel gain $\mu^{(i)}$ as well as the NPI variance $\nu^{2}_{i}$, we can write the $t^{\text{th}}$ symbol estimate of the $i^{\text{th}}$ user as follows:

Here, $\mathbf{W}^{(i,:)}=\left[\mathbf{W}^{(i,1)},\ldots,\mathbf{W}^{(i,{B})}\right]$ is a horizontal concatenation of the $i^{\text{th}}$ block row of (diagonal) submatrices of $\mathbf{W}$. The row vector $\mathbf{f}^{H}_{t}$ corresponds to the $t^{\text{th}}$ row of the IDFT matrix $\mathbf{F}_{L}^{H}$. Let $\mathbf{H}^{(:,j)}=\!\left[\mathbf{H}^{(1,j)},\ldots,\mathbf{H}^{({B},j)}\right]^{T}$ be the horizontal concatenation of the $j^{\text{th}}$ block column of (diagonal) submatrices of $\mathbf{H}$, consisting of the frequency-domain channel responses between the receive antennas and the transmit antenna associated with the $j^{\text{th}}$ user. We first compute the effective channel gain:

Since $\mathbf{W}^{(i,j)}$ and $\mathbf{H}^{(i,j)}$ are both diagonal matrices, we can write $\mu^{(i)}$ as a sum of per-subcarrier operations. In particular, let $\mathbf{w}^{H}_{i,w}$ be the $i^{\text{th}}$ row of $\mathbf{W}_{w}$ and $\mathbf{h}_{i,w}$ be the $i^{\text{th}}$ column of $\mathbf{H}_{w}$. Then, we obtain the effective channel gain as

We next compute the post-equalization NPI variance $\nu_{i}^{2}$ of the residual noise plus interference as

The MMSE equalization matrix can be written in two ways , i.e., either

Hence, we have $\mathbf{W}\!\left(E_{s}\mathbf{H}\mathbf{H}^{H}+N_{0}\mathbf{I}_{L{B}\times L{B}}\right)=E_{s}\mathbf{H}^{H}$; this allows us to rewrite the post-equalization NPI in compact form as follows :

We emphasize that both parameters $\mu^{(i)}$ and $\nu^{2}_{i}$ are functions of $\mathbf{A}^{-1}_{w}$, $\forall w$. Consequently, an explicit computation of the inverses $\mathbf{A}_{w}^{-1}$, $\forall w$, is necessary for the computation of LLR values using the approach detailed above.

III Approximate MMSE Detection via Neumann Series Expansion

The computation of all per-subcarrier inverses $\mathbf{A}^{-1}_{w}$, $\forall w$, in (1) is responsible for the main computational complexity of linear MMSE detection in SC-FDMA-based large-scale MIMO systems. For a conventional small-scale LTE uplink scenario, i.e., where the number of receive antennas ${B}$ and users ${U}$ is small (on the order of ${U},{B}\leq 6$), existing VLSI designs for linear detection, such as , compute the exact inverse explicitly. For large-scale MIMO systems with a large number of users ${U}$ however, the computation of the inverse $\mathbf{A}_{w}^{-1}$ can quickly result in excessive complexity. Hence, practical solutions for large-scale MIMO detection in LTE necessitate low-complexity matrix inversion methods—a corresponding approximate solution is proposed next.

For large-scale MIMO systems, where the number of receive antennas is larger than the number of single-antenna users, i.e., for ${U}\ll{B}$, the Gram matrices $\mathbf{G}_{w}$, and, consequently $\mathbf{A}_{w}$, become diagonally dominant . In fact, for i.i.d. Gaussian channel matrices $\mathbf{H}_{w}$ (with properly normalized entries) and in the large antenna limit, shows that $\mathbf{G}_{w}\rightarrow\mathbf{I}_{U}$. Inspired by this central property of large-scale MIMO, one can derive a low-complexity approximation of the inverse. In particular, let $\mathbf{A}_{w}\approx\mathbf{D}_{w}$, where $\mathbf{D}_{w}$ is the main diagonal of $\mathbf{A}_{w}$. As a result, the inverse $\mathbf{A}^{-1}_{w}$ can be approximated by $\mathbf{D}_{w}^{-1}$, which requires evidently much lower complexity than that of the exact inverse. Unfortunately, for realistic antenna/user configurations, such a crude approximation would cause a significant performance loss. Hence, to arrive at an accurate approximation of the inverse at low computational complexity, we propose to use a Neumann series expansion.

We start by rewriting the inverse $\mathbf{A}_{w}^{-1}$ with the following Neumann series expansion :

which holds if $\lim_{n\to\infty}(\mathbf{I}-\mathbf{X}^{-1}\mathbf{A}_{w})^{n}=\mathbf{0}_{{U}\times{U}}$ is satisfied. By decomposing the regularized Gram matrix $\mathbf{A}_{w}$ such that $\mathbf{A}_{w}=\mathbf{D}_{w}+\mathbf{E}_{w}$, where $\mathbf{D}_{w}$ is the main diagonal of $\mathbf{A}_{w}$ and $\mathbf{E}_{w}$ is the hollow, regularized Gram matrix, we can rewrite the Neumann series in (5) as

where we substitute $\mathbf{X}$ in (5) by $\mathbf{D}_{w}$. Note that if $\lim_{n\to\infty}(-\mathbf{D}^{-1}_{w}\mathbf{E}_{w})^{n}=\mathbf{0}_{{U}\times{U}}$, then the series expansion in (6) is guaranteed to converge.

The key idea of the proposed approximate inversion method is to keep only the first $K$ terms of the Neumann series (6). Concretely, we compute a $K$-term approximation as follows:

which can be computed at low computational complexity for approximations consisting of only a few Neumann series terms, i.e., for small values of $K$. With this approximation, the resulting approximate MMSE equalization matrix is given by $\widetilde{\mathbf{W}}_{w\,|\,K}=\widetilde{\mathbf{A}}^{-1}_{w\,|\,K}\mathbf{H}_{w}^{H}$. For $K=1$, we obtain $\widetilde{\mathbf{A}}^{-1}_{w\,|\,1}=\mathbf{D}_{w}^{-1}$, which is simply a scaled version of the MF detector, as $\widetilde{\mathbf{W}}^{-1}_{w\,|\,1}=\mathbf{D}_{w}^{-1}\mathbf{H}_{w}^{H}$. We emphasize that the row-wise scaling induced by $\mathbf{D}_{w}^{-1}$ does not affect the detection process, as long as $\mathbf{D}_{w}^{-1}$ exists. Hence, the proposed approximation (7) simply coincides with the MF detector for $K=1$. For $K=2$, we obtain $\widetilde{\mathbf{A}}^{-1}_{w\,|\,2}=\mathbf{D}_{w}^{-1}-\mathbf{D}_{w}^{-1}\mathbf{E}_{w}\mathbf{D}_{w}^{-1}$, whose computational complexity only scales with $O({U}^{2})$ operations; this is in contrast to the $O({U}^{3})$ complexity scaling required by computing an exact inverse. Hence, a second-order Neumann series approximation can be obtained at lower computational complexity. For $K=3$, we obtain

whose complexity scales with $O({U}^{3})$, which is equivalent to that of an exact inverse. Nevertheless, evaluating (8) requires fewer arithmetic operations than an explicit evaluation of $\mathbf{A}^{-1}$. Note that for $K\geq 4$, computing the exact inverse can be of lower complexity than the proposed approximation, e.g., when using a Cholesky factorization (see Section III-D).For approximations with $K=2^{n}$ terms and $n\geq 2$, efficient ways of evaluating (7) exist. In particular, a clever re-arrangement and factorization of terms yields solutions which only require $2(n-1)$ matrix multiplications.

III-B Analysis of the Approximation Error

We next analytically characterize the error induced by the approximate inverse (7) for MMSE estimation. To this end, we define the approximation error as $\mathbf{\Delta}_{w\,|\,K}=\mathbf{A}_{w}^{-1}-\widetilde{\mathbf{A}}_{w\,|\,K}^{-1}$, which is equivalent to

Now, consider the situation of using the approximate $\widetilde{\mathbf{A}}^{-1}_{w\,|\,K}$ in place of $\mathbf{A}_{w}^{-1}$ to compute the equalized frequency-domain symbols, i.e.,

is satisfied, then the approximation error approaches zero exponentially fast as $K\to\infty$. Moreover, one can show that (10) is a sufficient condition for (6) to converge.

We emphasize that this theorem provides conditionsThe result in (1) also holds for the case where the regularization term $N_{0}E_{s}^{-1}$ vanishes, which coincides to ZF detection. As a consequence, the condition (1) is rather pessimistic and is likely to be sub-optimal, especially for $N_{0}E_{s}^{-1}>0$. The derivation of a tighter condition is left for future work. for which the Neumann series converges with a certain probability; this can be accomplished by setting $\alpha=1$ and $K=1$ and by inspecting the convergence condition (10). Furthermore, Theorem 1 provides conditions for which the residual estimation error (9) is small. In both cases, we can see from Theorem 1 that increasing the ratio between the number of BS antennas ${B}$ and the number of users ${U}$ increases the probability of convergence. Moreover, for $\alpha<1$, increasing $K$ also increases the probability that the residual estimation error caused by a $K$-term approximation in (9) is smaller than $\alpha$.

III-C Channel Gain and NPI Variance Computation

Since $\mathbf{H}^{H}(E_{s}\mathbf{H}\mathbf{H}^{H}+N_{0}\mathbf{I}_{L{B}})=(E_{s}\mathbf{H}^{H}\mathbf{H}+N_{0}\mathbf{I}_{L{U}})\mathbf{H}^{H}$, we have:

Since (i) $a_{w}^{i,j}=g_{w}^{i,j}$, $\forall i\neq j$, (ii) $d^{(i,i)}_{w}=a^{(i,i)}_{w}$, and (iii) $d^{(i,i)}_{w}\gg a^{(i,j)}_{w}$ in the case where ${U}\ll{B}$, we can use the approximation $(d^{(i,i)}_{w})^{-2}\mathbf{a}_{i,w}^{H}\mathbf{g}_{i,w}\approx(d^{(i,i)}_{w})^{-1}g_{w}^{(i,i)}$. Hence, we propose the following low-complexity NPI approximation:

Note that our own simulations show that the low-complexity NPI approximation (15) performs well compared to the exact NPI variance (13). For example, the performance loss caused by the approximation compared to the exact NPI computation for ${U}=4$, ${B}=8$, and $K=3$ is less than $0.02$ dB at a BLER of $10^{-2}$ (cf. Section III-D2 for the simulation settings).

III-D Simulation Results

We next demonstrate the advantages and limitations of the proposed approximate matrix inversion approach in terms of computational complexity and error-rate performance. To assess the error-rate performance for practically relevant antenna configurations, we note that the Samsung Full-Dimensional MIMO prototype consists of 64 BS antennas, whereas the massive MIMO research platform developed at Rice University currently consists of 96 BS antennas (with plans for larger array sizes). Hence, we focus our results on the following cases: ${B}=64$, ${B}=128$, and ${B}=256$.

To demonstrate that the proposed approximate inverse exhibits (often significantly) lower complexity than an exact inverse, we chose a Cholesky decomposition-based inverse as a reference (see Section IV-E for algorithm details), as this method exhibits lower complexity compared to other inversion algorithms, including (but not limited to) direct matrix inversion, QR decomposition, or LU factorization . The computational complexity (characterized by the sum of real-valued divisionThe number of divisions is not significant for the total operation count., addition, and multiplicationTo obtain the real-valued multiplication count, we assumed four real-valued multiplications per one complex-valued multiplication. One could further reduce the number of the real-valued multiplications by using strength-reduction; this approach, however, maintains the trends observed in Fig. 1. operations) of an exact Cholesky-based inverse scales with $O({U}^{3})$, whereas the complexity of a $K=1$ and $K=2$ Neumann series expansion scales only with $O({U})$ and $O({U}^{2})$, respectively. The computational complexity of $K\geq 3$ is dominated by matrix-by-matrix multiplications, where the number of such operations grows linearly with $K$. For example, $K=3$ requires one matrix-by-matrix multiplication, whereas $K=4$ requires two. In general, a $K\geq 3$ term approximation requires $K-2$ matrix-by-matrix multiplications. As a result, the complexity of a $K\geq 3$ term approximation is $O((K-2)U^{3})$. Hence, we have $O({U}^{3})$ for $K=3$, which is equivalent to that of an exact Cholesky-based inverse. Consequently, a Neumann series approximation with $K\geq 3$ does not appear to be advantageous.

The overall operation counts of both methods are dominated by the number of real-valued multiplications and additions, where real-valued multiplication is more expensive than real-valued addition. Since asymptotic complexity scalings do not, in general, reveal the full truth, we count the number of real-valued multiplications of both methods in Figure 1 for varying numbers of users ${U}$. We observe that for $K\leq 3$, the Neumann series approach results in substantially lower complexity than the exact inversion approach. As expected, $K\geq 4$ results in higher complexity than a Cholesky-based exact inversion.

III-D2 Error-rate performance

Evidently, the reduction in complexity for $K\leq 3$ Neumann series terms comes at the cost of an approximation error (cf. Section III-B). To characterize the associated performance loss, we now compare the error-rate performance of the proposed approximate matrix inverse with the error-rate performance of the exact inversion for an LTE-based large-scale MIMO uplink system. To this end, we show simulation results of an SC-FDMA LTE uplink system with ${B}$ antennas at the BS and ${U}\leq{B}$ single-antenna users. In particular, we study a challenging communication scenario (from an error-rate perspective) and focus on the MCS (modulation and coding scheme) of the highest rate (i.e., MCS 28) and $20$ MHz bandwidth with $1200$ subcarriers, as specified by the LTE standard ; this mode corresponds to $64$-QAM, and a rate $\approx 0.75$ 3GPP LTE turbo code. In order to generate channel matrices that reflect a potentialTo the best of our knowledge, no specific channel model for large-scale MIMO systems is available in the open literature. real-world scenario, we use the WINNER-Phase-2 model . In addition, we assume a linear antenna array with an antenna spacing of $10/128\approx 0.0781$ m, which resembles that of the real-world channel measurement campaign in . At the BS, we use the exact and approximate soft-output MMSE detectors detailed above. Furthermore, we use a log-MAP LTE turbo decoder that performs $16$ (full-)iterations. Further, we define signal-to-noise-ratio (SNR) as ${B}E_{s}/N_{0}$, which corresponds to the average SNR per receive antenna.

Figures 2, 2, and 2 show the block-error rate (BLER) performance of the proposed approximate detection algorithm compared to that of an exact MMSE detector for ${U}=4$, ${U}=8$ and ${U}=12$, respectively.

We see that for small ratios between BS antennas and users, the MF detector (equivalent to $K=1$) and the Neumann series approximation for $K=2$ result in large residual errors.Compared to lower modulation orders, such as $16$-QAM (not shown here), $64$-QAM requires a relatively high SNR to perform well. Hence, considering the $10$% BLER requirement for LTE , the MF detector and $K=2$ term approximation are not suitable in practice in the considered $64$-QAM cases (note that this fact is also reflected by Theorem 1). For a larger number of BS antennas, this error floor can be recovered partially. Our own simulations have shown that the MF detector achieves $<10^{-2}$ BLER for ${U}=4$ and ${B}=512$. Furthermore, for $16$-QAM, our approximation method requires smaller values of $K$ (see for corresponding 16-QAM simulations in a large-scale MIMO-OFDM setting).

We see that for $64$-QAM, the proposed approximate inversion method with $K=3$ terms is able to approach the performance of the exact detector, i.e., the BLER performance loss is less than $0.25$ dB SNR at $10^{-2}$ BLER in all of the $K=3,{U}=4$ cases and the $K=3,{U}=8,{B}=256$ case. Hence, the proposed approximate inverse for $K=3$ can deliver the performance of an exact inversion at (often substantially) lower complexity for large ratios between BS antennas and users. For small antenna ratios, however, the approximate inverse with $K=3$ exhibits an error floor.

We conclude that systems with small ratios between BS and user antennas will need to resort to an exact inverse, while systems with large ratios can take advantage of the proposed approximate inverse. Hence, we next propose corresponding MIMO detection architectures for both, the approximate inverse and an exact Cholesky-based inverse.

IV VLSI Architecture

We now detail two VLSI architectures suitable for large-scale MIMO detection in 3GPP LTE-A. The first design implements the proposed approximate inversion approach and the second design implements an exact inverse; this enables us to perform a fair hardware complexity vs. error rate performance comparison (see Section V for the comparison).

The proposed general architecture is depicted in Figure 3 and consists of the following parts. The preprocessing unit performs matched filter computation, i.e., computes $\mathbf{y}^{\text{MF}}_{w}=\mathbf{H}_{w}^{H}\mathbf{y}_{w}$, the regularized Gram matrix, and the (approximate) inverse. Note that for the approximate inversion unit, we also output $\mathbf{D}_{w}^{-1}$ and $\mathbf{G}_{w}$, which are needed to compute the SINR (cf. Section III-C). To achieve the peak throughput specified in LTE-A , while being able to handle the (worst) case where the channel estimates change from subcarrier to subcarrier and from SC-FDMA symbol to SC-FDMA symbol (see, e.g, ), we use multiple instances of the preprocessing unit.In many practical scenarios, the channel estimates may change only slowly. Hence, one does not need to compute the inverse for every SC-FDMA symbol. This fact could be either exploited to reduce the power consumption or to increase the achievable throughput of our detector designs. The matched filter output, the (approximate) inverse, and the regularized Gram matrix, are then passed to the subcarrier processing unit. This unit performs equalization, i.e., computes $\hat{\mathbf{s}}_{w}=\mathbf{A}^{-1}_{w}\mathbf{y}^{\text{MF}}_{w}$ and the post-equalization SINR (detailed in Section II-B for the exact inverse and in Section III-C for the Neumann series approximation). To perform per-user data detection, a buffer is required that aggregates all equalized symbols and SINR values, which are computed on a per-subcarrier basis. The architecture then performs an IFFT, which transforms the equalized symbols from the subcarrier domain into the user domain (or time domain). The LLR computation unit finally computes, together with the buffered post-equalization NPI values, soft-output information in the form of max-log LLRs (2). We next provide the details for the key blocks of the proposed detector architecture.

IV-B Approximate Inversion and Matched Filter Units

In order to achieve high throughput, we propose a single systolic array that computes both, the regularized Gram matrix and the approximate inverse in four phases. The proposed architecture is detailed in Figure 4 and is capable of computing inverses for various $K$-term expansions, i.e., the number of Neumann series terms can be selected at run-time. As shown in Figure 4, the lower triangular systolic array consists of two distinct processing elements (PEs): (i) PEs on the main diagonal of the systolic array (referred to as PE-D) and PEs on the off-diagonal (referred to as PE-OD). As detailed next, both PEs have different modes in the four computation phases.

In the first phase, the ${U}\times{U}$ normalized regularized Gram matrix $\mathbf{A}_{w}/{B}=(\mathbf{G}_{w}+N_{0}E_{s}^{-1}\mathbf{I}_{U})/{B}$ is computed in ${B}$ clock cycles. Since $\mathbf{A}_{w}$ is diagonally dominant with diagonal entries close to ${B}$, i.e., the number of BS antennas, we reduce its dynamic range by computing a normalized version, whose entries on the main diagonal are close to $1$ by the ‘scale down’ unit shown in Figure 4; this trick mitigates dynamic-range issues, which are common for matrix inversion circuits implemented with fixed-point arithmetic. The systolic array also computes $\mathbf{D}_{w}^{-1}{B}$ from the diagonal entries of $\mathbf{A}_{w}/{B}$. These entries are computed in reciprocal units (denoted by ‘inv’ in Figure 4) residing in the PE-D units. The results $\mathbf{D}_{w}^{-1}{B}$ and $\mathbf{E}_{w}/{B}$ are then stored in register files distributed in the systolic array.

In the second phase, the systolic array computes $-\mathbf{D}_{w}^{-1}\mathbf{E}_{w}$, by using the matrices $\mathbf{D}_{w}^{-1}{B}$ and $\mathbf{E}_{w}/{B}$ computed in the first phase. Since the matrix $-\mathbf{D}_{w}^{-1}\mathbf{E}_{w}$ is not Hermitian, the systolic array computes the upper- and lower-triangular parts of $-\mathbf{D}_{w}^{-1}\mathbf{E}_{w}$ separately. As $\mathbf{D}_{w}^{-1}$ is a diagonal matrix, computation of $-\mathbf{D}_{w}^{-1}\mathbf{E}_{w}$ only requires a series of scalar multiplications (rather than a matrix multiplication).

In the third phase, the systolic array computes the $K=2$ term Neumann series approximation, i.e., $\widetilde{\mathbf{A}}^{-1}_{w\,|\,2}{B}=(\mathbf{D}_{w}^{-1}{B}-\mathbf{D}_{w}^{-1}\mathbf{E}_{w}\mathbf{D}_{w}^{-1}{B})$. To this end, it is important to realize that the matrix $\mathbf{D}_{w}^{-1}{B}-\mathbf{D}_{w}^{-1}\mathbf{E}_{w}\mathbf{D}_{w}^{-1}{B}$ is Hermitian, implying that only the lower triangular part needs to be computed. Furthermore, since $\mathbf{D}_{w}^{-1}{B}$ is diagonal, computation of $-\mathbf{D}_{w}^{-1}\mathbf{E}_{w}\mathbf{D}_{w}^{-1}{B}$ only requires entry-wise multiplications (instead of costly matrix multiplications). These scalar multiplications are carried out by loading $\mathbf{D}_{w}^{-1}{B}$ and $-\mathbf{E}_{w}\mathbf{D}_{w}^{-1}$ into all PEs and performing a scalar multiplication to compute $\mathbf{D}_{w}^{-1}\mathbf{E}_{w}\mathbf{D}_{w}^{-1}{B}$. Then, we add $\mathbf{D}_{w}^{-1}{B}$ to the result in the diagonal PEs. The result of this phase, i.e., $\mathbf{D}_{w}^{-1}{B}-\mathbf{D}_{w}^{-1}\mathbf{E}_{w}\mathbf{D}_{w}^{-1}{B}$, is stored in the distributed register files.

In the fourth phase, the $K$-term Neumann series approximation is computed with the results residing in the distributed register files. In particular, the systolic array first performs a matrix multiplication of $-\mathbf{D}_{w}^{-1}\mathbf{E}_{w}$ with $\widetilde{\mathbf{A}}^{-1}_{w\,|\,K-1}{B}$, and then adds $\mathbf{D}_{w}^{-1}{B}$ to the diagonal PE. The resulting $K$-term approximation $\widetilde{\mathbf{A}}^{-1}_{w\,|\,K}{B}$ is then stored in the register files. This phase can be repeated for a configurable number of iterations, which allows us to compute an arbitrary $K$-term approximation.

IV-B2 Matched filter computation

The matched filter (MF) unit consists of a linear array of ${U}$ PEs. Each PE is associated with one row of the Hermitian matrix $\mathbf{H}_{w}^{H}$, and contains a single multiply accumulate unit (MAC) and a scaling unit to normalize the result to $\mathbf{y}^{\text{MF}}_{w}/{B}$. The MF unit reads a new entry of $\mathbf{y}_{w}$ every clock cycle, and multiplies it with the corresponding entries in $\mathbf{H}_{w}^{H}$ in each PE and then, adds it the previous results; the final result is then normalized by $1/{B}$.

IV-C Equalization and SINR Computation Units

The equalization unit consists of a linear array of ${U}$ MAC units, and reads the normalized approximate inverse $\widetilde{\mathbf{A}}^{-1}_{w\,|\,K}{B}$ and the $\mathbf{y}^{\text{MF}}_{w}/{B}$ from the matched filter unit. For each clock cycle, this unit takes one column of $\widetilde{\mathbf{A}}^{-1}_{w\,|\,K}{B}$, multiplies it with one element from $\mathbf{y}^{\text{MF}}_{w}/{B}$, and adds the scaled column to the previous results. The unit outputs an equalized symbol $\hat{\mathbf{s}}_{w}$ every ${U}$ clock cycles.

IV-C2 SINR computation unit

IV-D IFFT and LLR Computation Units

In order to transform the per-subcarrier data into the user (or time) domain, we deploy a single Xilinx Discrete Fourier Transform IP LogiCORE unit (see for the specifications). This unit supports all forward and inverse DFT modes specified in 3GPP LTE , but we only make use of its IDFT capabilities. The IFFT unit reads and outputs data in a serial manner. For an IFFT transform size of $1200$ subcarriers, the core can process a new set of data every $3779$ clock cycles. This FFT unit achieves more than $317$ MHz on a Virtex-7 XC7VX980T FPGA and hence, achieves a throughput beyond $600$ Mb/s for $8$ users, $64$-QAM, and 20MHz bandwidth.

IV-D2 LLR computation unit

The LLR computation unit (LCU) generates max-log soft output values given the effective channel gains $\mu^{(i)}$ from the IFFT block and the post-equalization SINR values $\rho^{2}_{i}$ obtained from the SINR block. Since LTE specifies Gray mappings for all modulation schemes (BPSK, QPSK, 16-QAM, and 64-QAM), one can simplify the computation of the max-log LLR values in (2) by rewriting $L^{(i)}_{t}(b)={\rho^{2}_{i}}\lambda_{b}(\hat{x}^{(i)}_{t})$ and realizing that $\lambda_{b}(\cdot)$ is a piecewise linear function that depends on the bit index (see for the details). To this end, the LCU first scales the real and imaginary parts of the equalized time-domain symbol with the reciprocal of the effective channel gain $1/\mu^{(i)}$. Then, it evaluates the piecewise linear function $\lambda_{b}(\hat{x}^{(i)}_{t})$ and scales the result with the post-equalization SINR $\rho^{2}_{i}$. The resulting max-log LLR value is then delivered to the output of the unit. In order to minimize the circuit area, the proposed architecture evaluates each piecewise linear function with logical shifts and additions only. The reciprocals are computed with a lookup table that is stored in B-RAM units (see for architectural details). A single instance of the resulting LCU is able to processes one symbol every clock cycle, resulting in a peak throughput of $1.89$ Gb/s for 64-QAM at $317$ MHz.

IV-E Reference Cholesky-based Inversion Unit

In order to enable a fair performance/complexity assessment of the proposed approximate matrix inversion unit, we also implemented a reference unit that performs an exact matrix inversion. This unit simply replaces the approximate inverse unit detailed in Section IV-B. We next summarize the used Cholesky-based inversion algorithm and then, outline the corresponding VLSI architecture.

In the proposed exact inversion unit, we compute $\mathbf{A}^{-1}_{w}$ in three steps: (i) we form the regularized Gram matrix $\mathbf{A}_{w}=\mathbf{G}_{w}+N_{0}E_{s}^{-1}\mathbf{I}_{U}$; (ii) we perform a Cholesky decomposition according to $\mathbf{A}_{w}=\mathbf{L}_{w}\mathbf{L}_{w}^{H}$, where $\mathbf{L}_{w}$ is a lower-triangular matrix with real-values on the main diagonal ; (iii) we compute the inverse $\mathbf{A}_{w}^{-1}$ using an efficient forward/backward substitution procedure proposed in . Specifically, we first solve $\mathbf{L}_{w}\mathbf{u}_{i}=\mathbf{e}_{i}$ for $\mathbf{u}_{i}$, $i=1,\ldots,{U}$, where $\mathbf{e}_{i}$ is the $i^{\text{th}}$ unit vector, via forward substitution. We then solve $\mathbf{L}^{H}_{w}\mathbf{v}_{i}=\mathbf{u}_{i}$ for $\mathbf{v}_{i}$, $i=1,\ldots,{U}$, via back substitution, which leads to the desired inverse $\mathbf{A}^{-1}_{w}=[\,\mathbf{v}_{1}\cdots\mathbf{v}_{{U}}\,]$. Note that this approach avoids a costly matrix-by-matrix multiplication, which would be needed by directly computing $\mathbf{A}^{-1}_{w}=(\mathbf{L}_{w}^{H})^{-1}\mathbf{L}_{w}^{-1}$.

IV-E2 Cholesky decomposition architecture

The VLSI architecture for the Cholesky-based inverse differs from the one in Section IV-B. In particular, we deploy three separate units that compute (i) the regularized Gram matrix, (ii) the exact inverse using the above algorithm, and (iii) a forward/backward substitution unit to compute the inverse $\mathbf{A}_{w}^{-1}$. All units are detailed next and separated by pipeline stages.

The regularized Gram matrix is computed as a sum of outer products, i.e., as $\mathbf{G}_{w}=\sum_{i=1}^{{B}}\mathbf{r}_{i}\mathbf{r}^{H}_{i}$, where $\mathbf{r}_{i}$ designates the $i^{\text{th}}$ row of $\mathbf{H}_{w}$. Since the Gram matrix is symmetric, it can be computed efficiently with a triangular systolic array of multiply and accumulate units (MACs), similar to the array detailed in Section IV-B. The Gram computation unit reads one row of $\mathbf{H}_{w}$ at a time and is able to output a Gram matrix every ${B}^{\text{th}}$ clock cycle. To obtain the regularized Gram matrix $\mathbf{A}_{w}$, we add $N_{0}E_{s}^{-1}$ to the diagonal of $\mathbf{G}_{w}$ in the final clock cycle.

We then perform the Cholesky decomposition of $\mathbf{A}_{w}$ with a lower-triangular systolic array to obtain the lower-triangular matrix $\mathbf{L}_{w}$. The systolic array consists of two distinct processing elements (PEs): (i) the PEs on the main diagonal and (ii) the PEs on the off-diagonal. The data flow is similar to the linear systolic array (the “obvious case”) proposed in . The difference is that our design processes an incoming column of $\mathbf{A}_{w}$ with multiple PEs, whereas an incoming column is processed with a single PE in . As a result, our design is able to achieve the peak throughput requirements of LTE-A. In our design, the pipeline of one column of PEs is $16$ stages deep and streams out one column of $\mathbf{L}_{w}$ every clock cycle (after a latency of $16({U}-1)$ clock cycles). Consequently, the achieved throughput corresponds to one Cholesky decomposition every ${U}$ clock cycles.

IV-E3 Forward/backward-substitution architecture

The forward/backward substitution unit (FBSU) receives a lower-triangular matrix $\mathbf{L}_{w}$ as input, and computes $\mathbf{A}_{w}^{-1}=(\mathbf{L}_{w}^{H})^{-1}\mathbf{L}_{w}^{-1}$ as outlined in Section IV-E1. The FBSU consists of three major components: (i) a forward substitution unit (FSU), which solves for $\mathbf{L}_{w}\mathbf{u}_{i}=\mathbf{e}_{i}$, (ii) a backward substitution unit (BSU), which solves for $\mathbf{L}_{w}^{H}\mathbf{v}_{i}=\mathbf{u}_{i}$, and (iii) a Hermitian transpose unit, which computes $\mathbf{L}_{w}^{H}$. Since the computations for the FSU and the BSU are symmetric, we implement the forward substitution architecture and re-use it for the backward substitution, by reversing the order of the columns of the matrix $\mathbf{L}_{w}^{H}$ and vector $\mathbf{u}_{i}$ before reading them into the BSU. To simplify notation, we assume that the equation to be solved by the forward substitution corresponds to $\mathbf{Lx=b}$ for some $\mathbf{x}$ and $\mathbf{b}$. Since the forward substitution of solving the equation $\mathbf{L}\mathbf{x}_{i}=\mathbf{b}_{i}$ for each $\mathbf{b}_{i}$ $(i=1,\ldots,U)$ is independent, we use $U$ processor elements (PEs) to solve for all $\mathbf{x}_{i}$ in parallel. Each PE is implemented using a fully pipelined architecture, which consists of $U$ stages of computation logic. Each stage contains two multiplexers, a complex-valued multiplier, and a complex-valued subtraction. In each stage, either $\Delta_{i}={b}_{i}-\sum_{j}{L}_{i,j}{x}_{j}$ or $\Delta_{i}/{L}_{i,i}$ is computed according to the control signals. Therefore, for an input matrix $\mathbf{L}_{w}$ of dimension $U$, the FSU uses $U^{2}$ complex-valued multipliers; the entire FBSU utilizes $2U^{2}$ complex-valued multipliers. The matrix conjugate unit is implemented using multiplexers and $U$ FIFOs (realized by on-chip B-RAMs in the FPGA). The conjugate matrix $\mathbf{L}_{w}^{H}$ is also reordered based on the pattern of the input sequence of the BSU.

V Implementation Results and Trade-offs

The approximate detection engine for 3GPP-LTE and the exact Cholesky-based detector have been implemented on a Xilinx Virtex-7 XC7VX980T FPGA. The fixed-point parameters, FPGA implementation results, and the associated performance/complexity trade-offs are presented next.

In order to minimize the hardware complexity, fixed-point arithmetic is used in the entire design. The associated fixed-point parameters were determined via extensive simulations. In the following, the word-lengths refer to the real or imaginary part of a complex-valued number.

The channel matrices $\mathbf{H}_{w}$, the receive-vectors $\mathbf{y}_{w}$, and the noise variance $N_{0}E_{s}^{-1}$, are all quantized to $15$ bit. The word-length of the output of the Gram matrix and inversion unit are also set to $15$ bit; equivalently, the matched filter unit has $15$ bit at the input and output. For both matrix inversion circuits, all multiplications have been mapped onto Xilinx DSP48 slices. In order to achieve sufficient precision at minimum implementation complexity, the MAC registers within the DSP48 units are set to $22$ bit. The LUT in the reciprocal unit consists of $1024$ addresses with $12$ bit outputs. Hence, it can be implemented efficiently using a single block-RAM (B-RAM) available on the FPGA. The equalizer module uses a $15$ bit input and its output, which is stored in the data buffer, is quantized to $12$ bit. The buffer stores (complex-valued) data for $1200$ subcarriers and ${U}$ users. The SINR computation module has a $15$ bit input and $12$ bit output. The input and output of the IFFT unit are $12$ bit; the precision of the internal multipliers is set to $18$ bit. The inputs of the LLR computation are quantized to $12$ bit and the computed LLRs are represented by $8$ bit.

The resulting fixed-point performance is shown in Figure 2 (labeled by ‘FP’) for $64\times 4$ and $128\times 8$ systems. As it can be seen, the fixed-point implementation is virtually indistinguishable from the floating-point golden model. In particular, the implementation loss is less than $0.05$ dB SNR at 10% BLER.

V-B FPGA Implementation Results

Table I summarizes the key (post-place-and-route) implementation results of the proposed approximate and exact soft-output data detector for LTE-based massive MIMO wireless systems. We parameterized the architecture for ${U}$ and ${B}$ to explore the impact on the required FPGA resources and the corresponding throughput. The implementation results for antenna configurations of $128\times 8$ and $64\times 4$ are detailed in Table I. In order to support $75$ Mb/s data rate for each LTE-A user in 20 MHz bandwidth, we use multiple instances of the preprocessing unit. Specifically, we used $8$ and $5$ instances of approximate matrix inversion units for the $128\times 8$ and $64\times 4$ system, respectively. For the exact inverse, we used $6$ and $3$ regularized Gram matrix units for the $128\times 8$ and $64\times 4$ system, respectively. In addition, we used one Cholesky decomposition unit and one forward and backward substitution unit for both cases to meet the data rate requirements.

As shown in Table I, all designs are capable of running at $317$ MHz and the critical path is the routing between different blocks of the detector. For the $128\times 8$ and $64\times 4$ systems, the proposed units can achieve $603$ Mb/s and $301$ Mb/s, respectively. For the $64\times 4$ system, the design meets the $300$ Mb/s peak data rate requirement specified in LTE-A with $4$ users and 20MHz bandwidth. In addition, our design can scale beyond LTE-A specifications, i.e., the proposed designs can support up to $8$ users and still achieve a $75$ Mb/s per-user requirement.

In terms used resources on the Virtex-7 XC7VX980T FPGA, the approximate soft-output data detector is smaller than the Cholesky-based unit. There are notable saving in logic slices and DSP48 units. For $64\times 4$, $K=3$ uses 56% fewer LUT slices and 29% fewer DSP48 units compared to that of the Cholesky-based unit. For $128\times 8$, $K=3$ uses 19% fewer LUT slices and 26% fewer DSP48 units compared to that of the Cholesky-based unit. We emphasize that the savings in hardware resources become significantly larger as the number of users ${U}$ increases.

V-C Performance/Complexity Trade-off

Based on the simulated BLER results in Figure 2 and the associated FPGA implementation results, we are now ready to characterize the error-rate performance vs. hardware complexity trade-offs associated with the detector containing the proposed approximate matrix inversion and the Cholesky-based exact inversion. To this end, we show the associated hardware complexity against the minimum SNR required to achieve 10% BLER in Figure 5. Since both designs are dominated by multipliers, we define the hardware complexity as the number of multipliers required to achieve a $75$ Mb/s per-user throughput.

From Figure 5, we observe that the hardware complexity of the Cholesky-based detector is larger than that of the approximate inversion circuit for $K=3$ and $K=2$. In addition, for large ratios between the number of BS antennas to the number of users ${B}/{U}$, we clearly see that the SNR performance of the approximate inverse with $K=3$ and the exact inverse are very similar. For small ratios ${B}/{U}$, however, the performance difference between the approximate inverse and the exact inverse is rather large, which is reflected in the analysis shown Section III-B. Hence, the ratio ${B}/{U}$ determines whether an approximate or exact inversion is beneficial in a practical large-scale MIMO system. Note that for $128\times 8$ and $64\times 4$, the approximate inverse with $K=2$ is unable to achieve 10% BLER (cf. Figure 2). We note that when considering 16-QAM modulation (rather than 64-QAM modulation, as shown here), the approximate inversion for $K=2$ is capable of achieving similar performance as the exact inverse (see for corresponding simulation results).

V-D Related FPGA Designs for Linear Data Detection

A host of FPGA designs for linear data detection in conventional (small-scale) MIMO systems have been proposed in the literature . Unfortunately, all these designs differ in various ways. First, the corresponding architectures rely on different matrix inversion algorithms, such as the QR decomposition , Gram-Schmidt orthogonalization , LU decomposition, direct matrix inversion , divide-and-conquer methods . Second, all FPGA implementations do not generate soft outputs, with the exception of . Third, the designs were implemented on different FPGA types.

Since the soft-output detector implementations proposed in this paper are for large-scale MIMO systems having hundreds of BS antennas and none of the small-scale MIMO detector designs in was implemented on a Xilinx Virtex-7 FPGA, a fair comparison of our design with the above-mentioned implementations is difficult. Hence, we decided to resort to the comparison with our own reference circuit, i.e., the Cholesky-based inverse, as shown in Section V-C.

VI Conclusions

We have proposed a new soft-output data detector for large-scale (or massive) MIMO-based 3GPP LTE-Advanced (LTE-A) systems. The proposed solution is capable of performing high throughput detection in single-carrier frequency division multiple access (SC-FDMA)-based large-scale MIMO systems equipped with hundreds of antennas at the base station (BS). In order to achieve low computational complexity, we have proposed a new approximate linear detector relying on a Neumann series approximation of the matrix inverse. We have designed two reference VLSI architectures, one relying on the approximate inverse, the other on an exact Cholesky-based matrix inversion. Both architectures have been successfully implemented on a state-of-the-art Xilinx Virtex-7 FPGA, are suitable for systems equipped with $128$ BS antennas or fewer while serving up to $8$ users, and achieve more than $600$ Mb/s, exceeding the peak data rates specified in the 3GPP LTE-A uplink for 20 MHz bandwidth. Our FPGA implementation results reveal that for systems with a large ratio between the number of BS antennas and the number of users, the approximate matrix inversion is able to significantly reduce the hardware implementation complexity (compared to that of the exact inversion) with only a slight error-rate performance degradation. For systems with small ratios between the number of BS antennas and the number of users (as it is the case in, e.g., conventional, small-scale MIMO systems) one must resort to an exact inverse in order to avoid poor error-rate performance. This behavior is in accordance with the analytical results we have developed for the approximate matrix inverse. In summary, our FPGA implementation results demonstrate the practical feasibility of high-throughput data detection for 3GPP LTE-based large-scale MIMO systems. We finally note that a corresponding high-throughput ASIC design has recently been published in .

There are many avenues for future work. The development of detection algorithms that are able to perform iterative detection and decoding (as, e.g., in ) in large-scale MIMO systems is left for future work. Furthermore, the design of high-performance, near-optimal detection methods (e.g., based on the algorithms in ) that require low computational complexity for large-dimensional antenna configurations and for SC-FDMA is a challenging open research problem.

Appendix A Proof of Theorem 1

To prove Theorem 1, we need the following three Lemmata.

Let the scalars $x^{(k)}$ and $y^{(k)}$ for $k=1,\ldots,B$ be i.i.d. circularly symmetric complex Gaussian with unit variance. Then, $\mathsf{E}\!\left[\left|\sum_{k=1}^{{B}}x^{(k)}y^{(k)}\right|^{4}\right]=2{B}({B}+1)$.

The above steps can be summarized as follows. After expanding the quadratic expression, the non-zero terms can be written as $|x^{(k)}|^{4}|y^{(k)}|^{4}$ and $|x^{(k)}|^{2}|y^{(k)}|^{2}$, where $k=1,\ldots,{B}$. Then, there are ${B}$ terms of the form $|x^{(k)}|^{4}|y^{(k)}|^{4}$ and ${{B}\choose 2}$ of the form $|x^{(k)}|^{2}|y^{(k)}|^{2}$. The facts that $\mathsf{E}\left[|x^{(k)}|^{4}\right]=\mathsf{E}\left[|y^{(k)}|^{4}\right]=2$ and $\mathsf{E}\left[|x^{(k)}|^{2}\right]=\mathsf{E}\left[|y^{(k)}|^{2}\right]=1$ concludes the proof. ∎

Let ${B}>4$ and $x^{(k)}$, $k=1,\ldots,{B}$ be i.i.d. circularly symmetric complex Gaussian with unit variance and $g=\sum_{k=1}^{{B}}\,|\,x^{(k)}\,|\,^{2}$. Then,

We first rewrite $g$ as $2^{-1}\sum_{k=1}^{2{B}}\,|\,s^{(k)}\,|\,^{2}$ where $s^{(k)}$, $k=1,\ldots,2{B}$, are i.i.d. zero-mean real-valued Gaussian with unit variance. Then, $2g^{-1}$ is an inverse chi-square random variable with $2{B}$ degrees of freedom. The inverse chi-square distribution with $2{B}$ degrees of freedom $\chi(2{B})$ corresponds to an inverse-Gamma distribution with $2B$ degrees-of-freedom. The $4^{\text{th}}$ moment of this inverse chi-square distribution is given by $\frac{1}{16}({B}-1)({B}-2)({B}-3)({B}-4)$ and, hence, we obtain (16). ∎

The regularized Gram matrix corresponds to $\mathbf{A}_{w}=\mathbf{D}_{w}+\mathbf{E}_{w}=\mathbf{G}_{w}+N_{0}{E_{s}}^{-1}\mathbf{I}_{{U}\times{U}}$. Thus, each element on the $i^{\text{th}}$ row and $j^{\text{th}}$ column of $\mathbf{A}_{w}$, $a_{w}^{(i,j)}$ can be written as:

with $g_{w}^{(i,j)}$ corresponding to the $i^{\text{th}}$ row and $j^{\text{th}}$ column of the Gram matrix $\mathbf{G}_{w}$. We now have the following inequality:

which is obtained by omitting the non-negative regularization term $N_{0}E_{s}^{-1}$. By applying the Cauchy-Schwarz inequality, we can bound $\mathsf{E}\left[\|\mathbf{D}_{w}^{-1}\mathbf{E}_{w}\|^{2}_{F}\right]$ from above as

Application of Lemmata 2 and 3 to the first and second expected values, respectively, we obtain

We are now in position to prove Theorem 1. To this end, we start by using Markov’s inequality to obtain the following straightforward inequality:

References