Turbo EP-Based Equalization: A Filter-Type Implementation

1 Turbo EP-based Equalization: a Filter-Type Implementation Irene Santos, Juan Jos´ e Murillo-Fuentes, Eva Arias-de-Reyna, and Pablo M. Olmos Abstract—We propose a novel filter-type equalizer to improve the solution of the linear minimum-mean squared-error (LMMSE) turbo equalizer, with computational complexity constrained to be quadratic in the filter length. When high-order modulations and/or large memory channels are used the optimal BCJR equalizer is unavailable, due to its computational complexity. In this scenario, the filter-type LMMSE turbo equalization exhibits a good performance compared to other approximations. In this paper, we show that this solution can be significantly improved by using expectation propagation (EP) in the estimation of the a posteriori probabilities. First, it yields a more accurate estimation of the extrinsic distribution to be sent to the channel decoder. Second, compared to other solutions based on EP the computational complexity of the proposed solution is constrained to be quadratic in the length of the finite impulse response (FIR). In addition, we review previous EP-based turbo equalization implementations. Instead of considering default uniform priors we exploit the outputs of the decoder. Some simulation results are included to show that this new EP-based filter remarkably outperforms the turbo approach of previous versions of the EP algorithm and also improves the LMMSE solution, with and without turbo equalization. Index Terms—Expectation propagation (EP), linear MMSE, low-complexity, turbo equalization, ISI, filter-type equalizer. I. INTRODUCTION MANY digital communication systems need to transmit over channels that are affected by inter-symbol interference (ISI). The equalizer produces a probabilistic estimation of the transmitted data given the vector of observations [1]. Significant improvements are found when the previous estimation is given to a probabilistic channel decoder [2]. Equalization can be done in the frequency domain to avoid complexity problems associated with the inverse of covariance matrices [3]. In addition, feeding the equalizer back again with the output of the decoder, iteratively, yields a turbo-equalization scheme that significantly reduces the overall error rate [4]–[6]. The BCJR algorithm [7] performs optimal turbo equalization under the maximum a posteriori (MAP) criterion. It provides a posteriori probability (APP) estimations given some I. Santos, J.J. Murillo-Fuentes and E. Arias-de-Reyna are with the Dept. Teor´ ıa de la Se˜ nal y Comunicaciones, Escuela T. Superior de Ingenier´ ıa, Universidad de Sevilla, Camino de los Descubrimiento s/n, 41092 Sevilla, Spain. E-mail: {irenesantos,murillo,earias}@us.es P. M. Olmos is with the Dept. Teor´ ıa de la Se˜ nal y Comunicaciones, Universidad Carlos III de Madrid, Avda. de la Universidad 30, 28911, Legan´ es (Madrid), Spain. He is also with the Instituto de Investigaci´ on Sanitaria Gregorio Mara˜ n´ on (IiSGM). E-mail: [email protected]. Manuscript received September 7, 2017; revised January 10, 2018 and March 27, 2018; accepted April 25, 2018. This work was partially funded by Spanish government (Ministerio de Econom´ ıa y Competitividad TEC201678434-C3-{2-3}-R and Juan de la Cierva Grant No. IJCI-2014-19150) and by the European Union (FEDER). a priori information about the transmitted data. However, its complexity grows exponentially with the length of the channel and the constellation size, becoming intractable for few taps and/or multilevel constellations. In this situation, approximated BCJR solutions, such as [8]–[11], can be used. They are based on a search over a simplified trellis with only Me states, yielding a complexity which is linear in this number of states. However, the performance of these approaches is quite dependent on the channel realization and the order of the constellation used. In addition, these approximated BCJR solutions degrade rapidly if the number of survivor paths does not grow according to the total number of states. For these reasons, filter-based equalizers are preferred [12]. A quite extended filter type equalizer in the literature is based on the well-known linear minimum-mean squared-error (LMMSE) algorithm [5], [13], [14]. This LMMSE filter is an appealing alternative where the BCJR is computationally unfeasible due to its robust performance with linear complexity in the frame length, N, and quadratic dependence with the window length, W. From a Bayesian point of view, the LMMSE algorithm obtains a Gaussian extrinsic distribution by replacing the discrete prior distribution of the transmitted symbols with a Gaussian prior. A more accurate estimation for the extrinsic distribution can be obtained by replacing the prior distributions with approximations of the probability distribution. This can be done by means of the expectation propagation (EP) algorithm. The EP approach projects the approximated posterior distribution into the family of Gaussians by matching its moments iteratively with the ones of the true posterior. This algorithm has been already successfully applied to multiple-input multiple-output (MIMO) systems [15] and low-density parity-check (LDPC) channel decoding [16], [17], among others. It has been also applied to turbo equalization in a message passing approach as a way to incorporate into the BP algorithm the discrete information coming from the channel decoder [18], [19]. These message passing methods reduce to the LMMSE estimation if no turbo equalization is employed. A different approach is proposed in [20], [21] under the name of block EP (BEP) where, rather than applying EP after the channel decoder, it is used within the equalizer to better approximate the posterior, outperforming previous solutions. The computational complexity of previous EP-based equalizers is large for long frame lengths or memories of the channel. Due to its block implementation, the complexity of the BEP is quadratic in the frame length, becoming intractable for large frames [21]. To overcome this drawback, a smoothing EP (SEP) implementation is proposed in [22], but its complexity 2 is cubic with the memory of the channel. Furthermore, due to their iterative procedure, their computational load is roughly Stimes the one of the LMMSE counterparts, where Sis the number of iterations used in the EP algorithm, typically around 10 [15], [20], [21]. Besides, in both, BEP and SEP, uniform discrete priors are assumed for the constellation of the modulations when computing the EP approximations, even within the turbo equalization iterations, while the use of information from the decoder remains unexplored. The results developed in this paper focus on improving these previous EP-based equalizers [21], [22] both in computational complexity and performance. First, we improve the prior information used in the equalizer once the turbo procedure starts, forcing the true discrete prior to be non-uniform in contrast to the uniform priors used by previous EP-based approaches. As a result, we achieve a performance improvement. Second, the computational complexity of the EP algorithm is reduced to roughly a third part of that in [21], by optimizing the choice of EP parameters. Third, and most important, a new filter-type EP solution is designed. This solution is constrained to have linear complexity in the frame length and quadratic in the filter length, i.e., it is endowed with the same complexity order than the LMMSE filter. The novel EP-based filter proposed outperforms the LMMSE algorithm with a robust behavior to changes in the constellation size and the channel realization, as the BEP and SEP approaches do [21], [22]. In the experiments included, we show that the EP filter solution greatly improves the LMMSE solution with and without turbo equalization, specifically we have 2 dB gains for a BPSK, 5 dB for the 8-PSK and 6-13 dB for 16 and 64-QAM, respectively. In comparison with previous EP approaches, the EP filter matches their performance with BPSK constellations, and outperforms them with gains of 2 dBs for 8-PSK and 4-5 dBs for 16 and 64-QAM. We study the extrinsic information transfer (EXIT) charts [5], [23] of our proposal for a BPSK, where the EP-based filter achieves the same performance as the BEP. The scope of this paper encompasses time domain equalization. Frequency domain equalization has received a lot of attention as it usually achieves a complexity reduction for the block-wise processing [3], [14], [24], [25]. For this reason, derivation of a frequency domain counterpart for the proposed EP based turbo-equalizer remains as a future research line. Another promising research route is the application to MIMO with channels with memory [3], [26], [27]. The paper is organized as follows. We first describe in Section II the model of the communication system at hand. Section III is devoted to develop a new implementation of the EP-based equalizer considering non-uniform priors and studies the optimal values for the parameters. In Section IV, we review the formulation for the LMMSE filter in turbo equalization and describe the novel EP filter-type solution proposed. In Section VI, we include several simulations to compare both EP and LMMSE approaches. We end with conclusions. Through the paper, we denote the i-th entry of a vector uas ui, its complex conjugate as u˚and its Hermitian transpose as uH. We define δpuiqas the delta function that takes value one if ui“0and zero in other case. We use CNpu:µ,Σq to denote a normal distribution of a random proper complex vector uwith mean vector µand covariance matrix Σ. II. SYSTEM MODEL The model of the communication system is depicted in Fig. 1, including turbo equalization at the receiver. There are three main blocks: transmitter, channel and turbo receiver. A. Transmitter The information bit sequence, a“ ra1, ..., aKsJwhere aiP t0,1u, is encoded into the coded bit vector b“ rb1, ..., bVsJ with a code rate equal to R“K{V. After permuting the bits with an interleaver, the codeword c“ rc1, ..., cVsJ is partitioned into Nblocks of length Q“log2pMq,c“ rc1, ..., cNsJwhere ck“ rck,1, ..., ck,Qs, and modulated with a complex M-ary constellation Aof size |A| “ M. These modulated symbols, u“ ru1, ..., uNsJ, where each component uk“Rpukq ` jIpukq P A, are transmitted over the channel. Hereafter, transmitted symbol energy and energy per bit are denoted as Esand Eb, respectively. B. Channel The channel is completely specified by the CIR, i.e., h“ rh1, ..., hLsJ, where Lis the number of taps, and is corrupted with AWGN whose noise variance, σ2 w, is known. Each k-th entry of the complex received signal y“ ry1, ..., yN`L´1sJis given by yk“ L ÿ j“1 hjuk´j`1`wk“hJuk:k´L`1`wk,(1) where wk„CN`wk: 0, σ2 w˘and uk“0for kă1and kąN. C. Turbo receiver When no information is available from the channel decoder, the posterior probability of the transmitted symbol vector u given the whole vector of observations yyields ppu|yq 9 ppy|uqppuq(2) where, assuming equiprobable symbols, the prior would be given by ppuq “ 1 M N ź k“1ÿ uPA δpuk´uq.(3) This prior matches with the definition given in [21] but, as explained below, it is just valid before the turbo procedure. In a turbo architecture the equalizer and decoder iteratively exchange information for the same set of received symbols [5], [14]. Traditionally, this exchange of information is done in terms of extrinsic probabilities in order to improve convergence and avoid instabilities. The extrinsic information at the output of the equalizer (see Fig. 1), pEpuk|yq, is computed so as to meet the turbo principle [13]. These probabilities, pEpuk|yq, are approximated when the optimal solution is intractable. We will denote the approximation by qEpukq. 3 Channel Coder śModulator H(z) Equalizer Demap ϕp¨q ś´1Channel Decoder ś Map aibtctuk wk ykpEpuk|yqLEpctqLEpbtqpai LDpbtq LDpctqpDpukq Transmitter Channel Turbo receiver Fig. 1: System model. The extrinsic distributions are demapped, LEpck,j q “ log řukPA|ck,j “0pEpuk|yq řukPA|ck,j “1pEpuk|yq,(4) deinterleaved and given to the decoder as extrinsic loglikelihood ratios, LEpbtq. The channel decoder computes an estimation of the information bits, p a, along with the extrinsic LLRs on the coded bits, computed as LDpbt|LEpbqq “ log ppbt“0|LEpbqq ppbt“1|LEpbqq´LEpbtq.(5) These extrinsic LLRs are interleaved, mapped again and given to the equalizer as updated priors, pDpu|LEpbqq, which are computed as pDpuk|LEpbqq “ ÿ uPA δpuk´uq Q ź j“1 pDpck,j “ϕjpuq|LEpbqq, (6) with ϕjpuqdenoting the j-th bit associated to the demapping of symbol u. This process is repeated iteratively for a given maximum number of iterations, T, or until convergence. Note that in Fig. 1 we have included the computation of the extrinsic information within the equalization and channel decoding blocks. Note also that, once the turbo procedure starts, the prior in (2) is conditioned on the input at the channel decoder and the symbols are not equiprobable anymore. In this situation, the posterior distribution computed by the equalizer is given by ppu|yq 9 ppy|uq N ź k“1 pDpuk|LEpbqq,(7) The true posterior distribution in (2) and (7) has complexity proportional to ML. When this complexity becomes intractable, we will approximate it, denoting it as qpuq. In the following, we omit the dependence on the input at the decoder, LEpbq, to keep the notation uncluttered in the rest of the paper. It is also uncluttered in Fig. 1. III. NON-UNIFORM BEP TURBO EQUALIZER EP [28]–[32] is a technique in Bayesian machine learning that approximates a (non-exponential) distribution with an exponential distribution whose moments match the true ones. In this paper, we focus on computing a Gaussian approximation for the posterior in (7), which is clearly non Gaussian due to the product of discrete priors in (6). As introduced in [21], this is done by iteratively updating an approximation within the Gaussian exponential family by replacing the non Gaussian prior terms in (2) by a product of Gaussians1, i.e., qr`spuq 9 ppy|uq N ź k“1 ˜pr`s Dpukq “CN`y:Hu, σ2 wI˘N ź k“1 CN´uk:mr`s k, ηr`s k¯(8) The marginalization of the resulting approximated Gaussian posterior distribution for the k-th transmitted symbol and `-th EP iteration yields qr`spukq „ CN´uk:µr`s k, s2r`s k¯(9) where µr`s k“mr`s k`(10) `ηr`s khk H´σ2 wI`Hdiagpηr`sqHH¯´1 py´Hmr`sq, s2r`s k“ηr`s k´η2r`s khk H´σ2 wI`Hdiagpηr`sqHH¯´1 hk, (11) His the N`L´1ˆNchannel matrix given by H“ » — — — — — — — — — — – h10. . . 0 . . ........ . . hL ...0 0...h1 . . ........ . . 0. . . 0hL fi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi fl (12) and hkis the k-th column of H(see Appendix A for the demonstration). At this point it is interesting to remark that (9)- (11) are completely equivalent to equations (15)-(17) in [21]. Here we developed the values of the mean and variance for each symbol while in [21] they were computed in block form. The current description is simpler because we only include the elements of the covariance matrix that are used during the execution of the algorithm, excluding the non-diagonal elements. 1Note that in [21] we used an alternative expression for (8) (an exponential distribution with parameters γk“mk{ηkand Λk“1{ηk). 4 The mean and variance parameters in (8) are initialized with the statistics from the channel decoder as mr1s k“ÿ uPA u¨pDpuk“uq,(13) ηr1s k“ÿ uPA pu´mr1s kq˚pu´mr1s kq ¨ pDpuk“uq.(14) Then, they are updated in parallel and iteratively by matching the moments of the following distributions qr`s EpukqpDpukq moment matching ÐÑ qr`s EpukqCN´uk:mr``1s k, ηr``1s k¯ (15) where qr`s Epukqis an extrinsic marginal distribution computed as2 qr`s Epukq “ qr`spukq{˜pr`s Dpukq “ CN´uk:zr`s k, v2r`s k¯(16) where zr`s k“µr`s kηr`s k´mr`s ks2r`s k ηr`s k´s2r`s k ,(17) v2r`s k“s2r`s kηr`s k ηr`s k´s2r`s k .(18) Note that this equalizer differs from the one in [21] because we used different definitions for the true prior, pDpukq. In the current manuscript, we considered non-uniform and discrete priors, given by (6), during the moment matching procedure in the equalizer, while in [21] uniform priors as in (3) were considered by default even after the turbo procedure. To increase the accuracy of the algorithm, a damping procedure follows the moment matching in (15). We have defined an algorithm, described in Algorithm 1, called Moment Matching and Damping that runs these two procedures. A. The nuBEP algorithm Algorithm 2 contains a detailed description of the whole EP procedure, where Sis the number of EP iterations while Tis the number of turbo iterations. Note that the difference with the approach in [21] lies in the definition of the prior distribution used during the moment matching procedure. In [21], we use an uniform distribution (denoted with the indicator function), forcing the same a priori probability for the symbols, regardless of the information fed back from the decoder, during the moment matching employed in the equalizer even after the turbo procedure starts. In this paper, we refine the definition of the prior used in the moment matching of EP algorithm as in (6), considering non-uniform priors once the turbo procedure starts. For this reason, we named this algorithm non-uniform BEP (nuBEP) turbo equalizer. 2Note that in this paper we used qr`s Epukqto denote the extrinsic marginal distribution, while in [21] we denoted as qr`szkpukqand called it cavity marginal function. Algorithm 1 Moment Matching and Damping Given inputs:µr`s pk, σ2r`s pk, zr`s k, v2r`s k, mr`s k, ηr`s k 1) Run moment matching: Set the mean and variance of the unnormalized Gaussian distribution qr`s Epukq ¨ CN´uk:mr``1s k,new, ηr``1s k,new¯(19) equal to µr`s pkand σ2r`s pk, to get the solution ηr``1s k,new “σ2r`s pkv2r`s k v2r`s k´σ2r`s pk ,(20) mr``1s k,new “ηr``1s k,new˜µr`s pk σ2r`s pk ´zr`s k v2r`s k¸.(21) 2) Run damping: Update the values as ηr``1s k“˜β1 ηr``1s k,new ` p1´βq1 ηr`s k¸´1 ,(22) mr``1s k“ηr``1s k˜βmr``1s k,new ηr``1s k,new ` p1´βqmr`s k ηr`s k¸.(23) if ηr``1s kă0then ηr``1s k“ηr`s k, mr``1s k“mr`s k.(24) end if Output:ηr``1s k, mr``1s k B. On the election of EP parameters The moment matching condition explained in (15) determines the optimal operation point found by the EP approximation. By repeating this procedure, we allow to find a stationary solution for the operation point. In order to avoid instabilities and control the accuracy and speed of convergence, some EP parameters are introduced. These parameters are the number of EP iterations (S), a minimum allowed variance () and a damping factor (β). Based on recent studies, these EP parameters can be further optimized [33]–[35]. Following the guidelines in those papers and after extensive experimentation, in the general case we found out that instabilities can be controlled by setting3“1e´8. Regarding the accuracy of the algorithm, it is convenient to start with a conservative value of the damping parameter βin Algorithm 2. The value β“0.1forces our algorithm to move slowly towards the EP solution. Once the turbo procedure starts, we let the damping parameter grow in order to speed up the achievement of the EP solution, reducing the value of Sfrom 10 in [21] to 3. A simple rule for determination of βthat fulfills this requirements and leads to good performance is an exponential growth with a saturation value of 0.7, i.e., β“minpexpt{1.5{10q,0.7q, where tP r0, Tsis the number of the current turbo iteration. With this criterion the number of EP iterations after the turbo procedure starts is reduced to S“3, hence reducing the computational complexity by more than a third. 3Parameters have been chosen to optimize turbo equalization [35]. 5 Algorithm 2 nuBEP Turbo Equalizer Initialization: Set pDpukq “ 1 MřuPAδpuk´uqfor k“ 1,...,N for t“1, ..., T do 1) Compute the mean mr1s kand variance ηr1s kgiven by (13) and (14), respectively. for `“1, ..., S do for k“1, ..., Ndo 2) Compute the k-th extrinsic distribution as in (16), i.e., qr`s Epukq “ CN´uk:zr`s k, v2r`s k¯(25) where zr`s kand v2r`s kare given by (17) and (18), respectively. 3) Obtain the distribution p pr`spukq 9 qr`s EpukqpDpukqand estimate its mean µr`s pkand variance σ2r`s pk. Set a minimum allowed variance as σ2r`s pk“maxp, σ2r`s pkq. 4) Run the moment matching and damping procedures by executing Algorithm 1. end for end for 5) With the values mrS`1s k, ηrS`1s kobtained after the EP algorithm, calculate the extrinsic distribution qEpukq. 6) Demap the extrinsic distribution and compute the extrinsic LLR, LEpck,jq, by means of (4). 7) Run the channel decoder to output pDpukq end for Output: Deliver LEpck,jqto the channel decoder for k“ 1,...,Nand j“1, . . . , Q IV. FILTER-TYPE TURBO EQUALIZATION A. LMMSE filter In this subsection we review the formulation of the LMMSE-based filter [5], [13], [14], modified to allow for unnormalized transmitted energy and a different computation of the extrinsic distribution. The LMMSE-based filter [5], [13], [14] estimates one symbol per k-th iteration, uk, given aW-size window of observations, yk“ ryk´W2, ..., yk`W1sJ, where W“W1`W2`1. This procedure differs from [36], where each transmitted symbol is estimated given the whole vector of observations, y. The LMMSE equalizer approximates the prior for each symbol, pDpukq, as a Gaussian pDpukq « ˜pDpukq “ CNpuk:mk, ηkq,(26) where the mean, mk, and variance, vk, are a priori statistics for each transmitted symbol, given by (13) and (14), respectively. For the first iteration of the turbo equalization no a priori information is available and a suitable initialization is mk“0, ηk“Es, which boils down to mk“0,ηk“1when normalizing the energy [5], [13], [14]. Given the current prior and the channel impulse response (CIR), the LMMSE filter computes a Gaussian approximation of the posterior probability of each symbol. When a turbo scheme is used, the equalizer and decoder exchange extrinsic information [6]. Through the turbo equalization iterations, the a priori statistics in (26) are updated with the information fed back from the channel decoder. Rather than computing the posterior distribution as in (7), the LMMSE filter [5] considers the a posteriori probabilities with respect to the estimated transmitted symbol, puk. For this reason, and to keep the same notation than in [5], we will denote the approximated posterior as qpuk|pukq. With this posterior distribution in mind, the extrinsic probability at the output of the LMMSE filter can be computed as qEpuk|pukq “ qpuk|pukq ˜pDpukq.(27) This distribution is Gaussian and can be derived from the extrinsic distribution of the estimated symbol computed in [5], as shown in Appendix B, yielding qEpuk|pukq “ CN`uk:zk, v2 k˘(28) where zk“ck Hpyk´HWmk`mkhWq ckHhW ,(29) v2 k“ck HhWEsp1´hW Hckq pckHhWq2,(30) and, in turn, ck“´Σk` pEs´ηkqhWhW H¯´1 EshW,(31) HW“» — — — — – hL. . . h10 ...... ...... 0hL. . . h1 fi ffi ffi ffi ffi fl (32) is the Wˆ pW`L´1qchannel matrix, hWis the (W2`L)-th column of HWand mk“ rmk´L´W2`1, ..., mk`W1sJ,(33) Vk“diagpηk´L´W2`1, ..., ηk`W1q,(34) Σk“σ2 wI`HWVkHW H.(35) The computational complexity is dominated by (31), which has to be recomputed every k-th iteration. Hence, the complexity is OpNW2q. This complexity can be further reduced by relying on some approximations proposed in [5], [14]. B. EP filter (EP-F) A novel EP filter-type is developed in this subsection to improve the accuracy and performance of the LMMSE-based filter explained above. As explained in Subsection IV-A, if the LMMSE filter is run, the prior of each symbol is approximated by a Gaussian with the statistics given by the decoder, i.e., with mean and variance given by (13) and (14), respectively. By using the EP algorithm we approximate the posterior distribution with a Gaussian family. Since the posterior distribution includes the true discrete priors, we take into account the discrete nature of symbols. At every iteration of the EP algorithm, `, we approximate the product of priors of individual symbols in (7) as a product 6 Algorithm 3 EP-F Initialization: Set pDpukq “ 1 MřuPAδpuk´uqfor k“ 1,...,N for t“1, ..., T do 1) Compute the mean mr1s kand variance ηr1s kgiven by (13) and (14), respectively. for `“1, ..., S do for k“1, ..., Ndo 2) Compute the k-th extrinsic distribution as in (28), i.e., qr`s Epuk|pukq “ CN´uk:zr`s k, v2r`s k¯(37) where zr`s kand v2r`s kare given by (29) and (30), respectively. 3) Obtain the distribution ppr`spukq 9 qr`s Epuk|pukqpDpukqand estimate its mean µr`s pkand variance σ2r`s pk. Set a minimum allowed variance as σ2r`s pk“maxp, σ2r`s pkq. 4) Run the moment matching and damping procedures by executing Algorithm 1. end for end for 5) With the values mrS`1s k, ηrS`1s kobtained after the EP algorithm, calculate the extrinsic distribution qEpuk|pukq in (28). 6) Demap the extrinsic distribution and compute the extrinsic LLR, LEpck,jq, by means of (4). 7) Run the channel decoder to output pDpukq end for Output: Deliver LEpck,jqto the channel decoder for k“ 1,...,Nand j“1, . . . , Q of NGaussians, ˜pr`s Dpukq “ CN´uk:mr`s k, ηr`s k¯, whose parameters (means and variances) are adjusted to find a better approximation, qr`spuq 9 ppy|uqśN k“1˜pr`s Dpukq, to the true posterior. Similarly to (27)-(28), for each k-th symbol, we first compute the current extrinsic distribution, qr`s Epuk|pukq “ qr`spuk|pukq ˜pr`s Dpukq.(36) Now, a more accurate posterior distribution can be obtained by finding a new Gaussian approximation, ˜pr``1s Dpukq, to match the moments of qr`s Epuk|p ukq˜pr``1s Dpukqand qr`s Epuk|p ukqpDpukq, as in (15). With these new values for the mean, mr``1s k, and variance, ηr``1s k, we can recompute a new extrinsic distribution qr``1s Epuk|pukq, which is more accurate than the one in (28). The final extrinsic distribution delivered to the decoder is the one obtained after the last iteration of the EP algorithm, following (36). We denote this new algorithm as EP-filter (EP-F). Algorithm 3 is a detailed description of its implementation. Note that the main difference between Algorithm 2 and Algorithm 3 lies in the computation of the extrinsic distribution, i.e., equations (25) and (37). The computational complexity is also dominated by (31), which has to be computed for each symbol and each `-th iteration. Hence, the complexity is Stimes the LMMSE complexity, i.e. OpSNW2q, where Sis the number of iterations of the EP-F. At this point, it is interesting to remark that the approximations proposed in [5], [14] to further reduce the complexity cannot be applied when the EP is used. The reason is that these approximations remove (at some points) the prior variance computed by the decoder, setting it to one. V. RELATION TO PREVIOUS APPROACHES A. Update of the priors We improve the prior information used in the equalizer once the turbo procedure starts, forcing the true discrete prior to be non-uniform in contrast to the uniform priors used by previous EP-based approaches. In previous proposals [20]–[22], the probabilities from the channel decoder, pDpukq, were used to initialize, at the beginning of every iteration of the turbo-equalization, the product of Gaussians that in the EP approximation replaces the product of priors, ˜pr1s Dpukq. But when the moment matching was performed in the EP algorithm, i.e., qr`s EpukqIukPA moment matching ÐÑ qr`s Epukq˜pr``1s Dpukq,(38) the true priors used were uniformly distributed following IukPA“1 Mÿ uPA δpuk´uq.(39) In the current proposal, we keep the initialization of the Gaussians in every step of the turbo-equalization, ˜pr1s Dpukq, but also propose to replace the uniform priors in (39) by nonuniform ones in the moment matching step, as explained in (15), i.e., qr`s EpukqpDpukq moment matching ÐÑ qr`s Epukq˜pr``1s Dpukq,(40) where pDpukq “ ÿ uPA δpuk´uq Q ź j“1 pDpck,j “ϕjpuqq,(41) Note that the different definition of priors -(39) in previous proposals, (41) in this manuscriptis the difference between the currently proposed nuBEP algorithm and the BEP in [21], with remarkable improvements. B. Parameter Optimization The computational complexity of the EP algorithm is reduced to roughly a third part of that in [21], by optimizing the choice of EP parameters. In particular, we propose some new values for and β, that control numerical instabilities in the EP updates, and S, the number of iterations of the EP equalizer. The parameters proposed in this paper reduce the number of iterations in turbo equalization to S“3, rather than the S“10 iterations that were used in [21]. 7 C. Filter-type solution The new filter-type EP solution proposed is constrained to have linear complexity in the frame length and quadratic in the filter length, i.e., it is endowed with the same complexity order than the LMMSE filter. This complexity is not quadratic with the block length as the one of the BEP [21] nor cubic with the window length as complexity of the SEP [22]. D. Equalization solved with EP Regarding the EP-based equalizers proposed by other authors, the approach in [18], [19] should be mentioned. These proposals deal just with how to pass information between the channel decoder and the LMMSE equalizer. Our proposal first focuses on the EP based equalization, performed independently of the turbo iterations. Therefore the approaches are quite different. Issues such as how to use the priors in the moment matching within the EP equalizer or the damping do not arise in these proposals where the improvement is related only to the handling of probabilities between blocks. VI. SIMULATION RESULTS In this section, we compare the performance of both the block LMMSE and EP-F equalizers for different scenarios. We also include the performance of the BEP [21] and the AWGN bound as references. Note that the MMSE filter [5] has not been included in the simulations since the block LMMSE exhibits equal or better performance than any filtering approaches based on the LMMSE algorithm. We did not include the SEP algorithm since it exhibits the same performance as the block implementation, as shown in [22]. We also include the nuBEP approach to illustrate the quite improved behavior when using non-uniform priors at each EP iteration, even reducing from 10 to 3 the number of iterations of the EP approach. The EP parameters have been selected as explained in Subsection III-B, both for the nuBEP and EPF methods. For a full performance comparison with BCJR approximations, such as M-BCJR [8], M*-BCJR [10], RSBCJR [9], NZ and NZ-OS [11], please see [21]. In Table I we include a detailed comparison of the complexity of all the simulated algorithms. Above we include the computational complexity of previous algorithms in [21] (BEP) and [22] (SEP), the block and filter implementation of the LMMSE and BCJR approaches. Below we provide the complexity for the new approaches in this paper, i.e., the proposed nuBEP and EP-F. Parameter Wis typically around two times the length of the channel, L. Here, we simulate the scenarios in [13], [14], using the same channel responses and modulations. Other modulations are also considered. The absolute value of LLRs given to the decoder is limited to 5 in order to avoid very confident probabilities. We use a (3,6)-regular LDPC of rate 1/2, and belief propagation as decoder with a maximum of 100 iterations. The window length in the filtered approach is set to W“W1`W2`1, where W1“2Land W2“L`1as suggested in [14]. In the following, we first include a section to analyze the performance of our approach in a low complexity scenario with BPSK modulation, similarly to [14]. The optimal BCJR TABLE I: Complexity comparison between algorithms. Algorithm Complexity per turbo iteration BCJR NML BEP 10LN2 block-LMMSE LN2 SEP 10NW3 LMMSE filter NW2 nuBEP 3LN2 EP-F 3NW2 algorithm can be run in this scenario with a low enough computational complexity and is used as bound. Next, we include a section to analyze the behavior of the algorithms in a large complexity scenario, where we use high-order modulations such as 8-PSK, 16-QAM and 64-QAM. A. BPSK scenario In Fig. 2 we include the BER, averaged over 104random frames, for the LMMSE, BEP [21], nuBEP, EP-F and BCJR equalizers with a BPSK modulation and two different channel responses and lengths of encoded words: h“ r0.227 0.46 0.688 0.46 0.227sJand V“4096 bits in Fig. 2 (a)-(c) and h“ r0.407 0.815 0.407sJand V“1024 bits in Fig. 2 (d)-(f). The channel responses were selected following the simulations in [13], [14]. The performances of blockalgorithms, BEP and nuBEP, are very similar to the equivalent forward filtering approach. When the nuBEP algorithm is applied, 2 and 1.5 dBs gains are obtained compared to LMMSE approach in the turbo scenario, for the two simulated scenarios, respectively. The EP-F exhibits a performance similar to that of the nuBEP. In Fig. 3 we include the EXIT charts of the BEP [21], nuBEP, EP-F, LMMSE and BCJR for the channel response h“ r0.227 0.46 0.688 0.46 0.227sJas in [5], [13], BPSK modulation with Eb{N0“9(solid) and 7dB (dashed). The EXIT chart of the LDPC encoder of 2048{4096 and R“1{2used is also depicted (solid). The horizontal and vertical axis depict the mutual information at the input, Ii, and the output, Io, respectively. We use arrows to show the evolution of the mutual information along the turbo iterations for Eb{No“9dB. Vertical (horizontal) arrows indicate the improvement in the mutual information each time the equalizer (channel decoder) is executed. When no a priori information is given to the decoder, i.e., Ii“0, both BEP and EP-F provide a higher value for the mutual information at the output, Io, than the LMMSE approach, i.e., they start from a more accurate estimation even before the turbo equalization. This greatly improves the performance as it enlarges the gap between the equalizer and the channel decoder EXIT curves. It can be seen that the LMMSE approach will fail when Eb{No“7dB, because both curves intersect. Note that the wide EXIT tunnel from the equalizer to the LDPC decoder is suggesting that the code is not optimum in terms of capacity [37]. An optimal code in this sense would 8 0 5 10 15 10´5 10´4 10´3 10´2 10´1 Eb{N0(dB) BER Channel h“ r0.227 0.46 0.688 0.46 0.227sJ (a) Without turbo feedback 0 5 10 15 10´5 10´4 10´3 10´2 10´1 Eb{N0(dB) BER (b) After two turbo loops 0 5 10 15 10´5 10´4 10´3 10´2 10´1 Eb{N0(dB) BER (c) After five turbo loops 2 4 6 8 10´5 10´4 10´3 10´2 10´1 Eb{N0(dB) BER Channel h“ r0.407 0.815 0.407sJ (d) Without turbo feedback 2 4 6 8 10´5 10´4 10´3 10´2 10´1 Eb{N0(dB) BER (e) After two turbo loops 2 4 6 8 10´5 10´4 10´3 10´2 10´1 Eb{N0(dB) BER (f) After five turbo loops Fig. 2: BER along Eb{N0for BEP [21] (˝), nuBEP (˝), EP-F (˚), block-LMMSE (Ź) and BCJR (˛) turbo equalizers, BPSK, codewords of V“4096 bits (a)-(c) and V“1024 bits (d)-(f) and two different channel responses. Black lines represent the AWGN bound. exhibit an EXIT chart near below the one of the equalizer and above the curve of the 1{2rate LDPC code used. The design of this code for the channel equalization response is out of the scope of this paper and remains as a future line of research. B. Large complexity scenario In Fig. 4 we simulate the same scenario of Fig. 2, but using an 8-PSK modulation rather than a BPSK. It can be observed that after increasing the order of the modulation, the EP-F approach presented in this work performs identically as its block counterpart, greatly improving the performance of the LMMSE algorithm before and after the turbo procedure. An improvement of the BER of the EP-F with respect to the one of the BEP approach in [21], after turbo equalization, can also be observed. In Fig. 5 we depict the BER performance after five turbo loops for channels h“ r0.227 0.46 0.688 0.46 0.227sJ in (a) and h“ r0.407 0.815 0.407sJin (b) with different modulations. We use solid lines to represent a 64-QAM constellation and dashed lines for a 16-QAM. We sent codewords of length V“4096 in both scenarios. It can be observed that the performance of the EP-F matches with the one of its block implementation proposed in this paper (nuBEP) when a 16-QAM is used. However, the EP-F approach slightly degrades with a 64-QAM, where the block nuBEP gets the most accurate performance. Note that the behavior of the EP-F could be improved by increasing the length of the filter, 9 0 0.2 0.4 0.6 0.8 1 0 0.2 0.4 0.6 0.8 1 Ii Io BEP [21] nuBEP EP-F LMMSE BCJR LDPC Decoder Fig. 3: EXIT charts for the decoder, the BEP [21], nuBEP, EP-F, LMMSE and BCJR equalizers with 7(dashed) and 9(solid) dB of Eb{N0, BPSK modulation and codewords of V“4096. yielding the performance of its block implementation. We have a remarkable improvement of 3-5dB with respect to the BEP in [21] and of 7-13 dB compared to the LMMSE algorithm. For the sake of completeness, we include Fig. 6 to show how the BER changes along the turbo iterations and different block lengths at Eb{N0“13 dB for an 8-PSK and h“ r0.227 0.46 0.688 0.46 0.227sJ. The nuBEP algorithm is represented in (a) and the filter approach EP-F in (b). It can be observed that BER is higher for shorter codes and it improves when the code length is increased, as expected. Also, the BER does not significantly improve after the fifth turbo iteration. In the view of these results we stopped running our turbo equalizers after five turbo iterations in the experiments presented above. VII. CONCLUSION In a previous work, we presented a novel equalizer based on expectation propagation (EP) [21]. This solution presents quite an improved performance compared to previous approaches in the literature, both for hard, soft and turbo detection. The solution was presented as a block-wise solution and it was therefore denoted as block-EP (BEP). The major advantage of the BEP lies in the fact that its computational complexity does not grow exponentially with the constellation size and channel memory, as opposed to most equalizers, which are unfeasible for moderate values of these parameters. However, it exhibits a quadratic increase with the size of the transmitted word, V. To avoid this problem, filter-type equalizers are usually preferred [12]. For this reason, we proposed a smoothing EP (SEP) equalizer in [22]. However, the SEP has a computational complexity cubic in the channel length, L. Both BEP and SEP equalizers make use of a moderate feedback in the sense that an initial uniform discrete prior is assumed at the beginning of each execution of the EP algorithm, even after the turbo procedure has started. In this paper, we first propose a design to include the non-uniform discrete nature of the priors from the decoder in the EP algorithm, which amounts to a stronger feedback, quite outperforming the previous BEP and SEP approaches. Second, we develop a reduced-complexity approach by proposing better values of the EP parameters. The resulting algorithm has been denoted as nuBEP, and it significantly outperforms the BEP reducing the computational complexity to less than the third part. Finally, we adapt the EP block equalizer to the filter-type form, emulating the Wiener MMSE filter-type [14]. Therefore, we mimic the structure of the filter-type MMSE equalizer. The EP is used to better approximate the posteriors of a windowed version of the inputs, shifted for every new output estimate. As a result, we present a novel solution dealing with Winputs at a time and with quadratic computational complexity in W. This novel solution, the EP-F, despite the reduction in the computational complexity, exhibits a performance in terms of BER quite close to that of its block counterpart, the nuBEP. Furthermore, it remarkably improves the performance of the LMMSE turboequalizer, with same complexity order in terms of Land V. In the included experiments, for channels usually used as benchmarks in the literature, gains in the range 5-13 dB are reported for 8-PSK, 16-QAM and 64-QAM modulations. One of the main benefits of this new proposal is to reduce the computational complexity, reducing it to be of quadratic order with the filter length. Other approaches, such as those solutions working on the frequency domain [3], could be investigated to achieve this goal. In this paper we face the equalization in single-input single-output channels, the application to MIMO channels with memory [27] remains unexplored. APPENDIX A PROOF OF (10) AND (11) In [21], the posterior distribution used for BEP is qr`spuq „ CN´u:µr`s,Σr`s¯(42) where µr`s“Σr`spσ´2 wHHy`diagpηr`sq´1mr`sq,(43) Σr`s“´σ´2 wHHH`diagpηr`sq´1¯´1 .(44) By a direct application of the Woodbury identity, equation (44) can be rewritten as Σr`s“diagpηr`sq ´ diagpηr`sqHHC´1Hdiagpηr`sq(45) where C“Hdiagpηr`sqHH`σ2 wI.(46) The k-th diagonal element of (45) yields (11). Regarding (43), it can be divided into two terms µr`s“Σr`sσ´2 wHHy loooooomoooooon T1 `Σr`sdiagpηr`sq´1mr`s loooooooooooomoooooooooooon T2 .(47) We apply the following identity [31], pA´1`BHD´1Bq´1BHD´1“ABHpBABH`Dq´1(48) to the first term, T1, in (47), yielding T1“diagpηr`sqHHC´1y.(49)