Maximum-Likelihood Semiblind Equalization of Doubly Selective Channels Using the EM Algorithm

Maximum-likelihood semi-blind joint channel estimation and equalization for doubly selective channels and single-carrier systems is proposed. We model the doubly selective channel as an FIR ﬁlter where each ﬁlter tap is modeled as a linear combination of basis functions. This channel description is then integrated in an iterative scheme based on the expectation-maximization (EM) principle that converges to the channel description vector estimation. We discuss the selection of the basis functions and compare various functions sets. To alleviate the problem of convergence to a local maximum, we propose an initialization scheme to the EM iterations based on a small number of pilot symbols. We further derive a pilot positioning scheme targeted to reduce the probability of convergence to a local maximum. Our pilot positioning analysis reveals that for high Doppler rates it is better to spread the pilots evenly throughout the data block (and not to group them) even for frequency-selective channels. The resulting equalization algorithm is shown to be superior over previously proposed equalization schemes and to perform in many cases close to the maximum-likelihood equalizer with perfect channel knowledge. Our proposed method is also suitable for coded systems and as a building block for Turbo equalization algorithms.


Introduction
Next generation cellular communication systems are required to support high data rate transmissions for highly mobile users.These requirements may lead to doubly selective channels, that is, channels that experience both frequency-selective fading and time-selective fading.The frequency selectivity of the channel stems from the requirement to support higher data rates that necessitates the usage of larger bandwidth.Time selectivity arises because of the need to support users traveling at high velocities as well as the usage of higher carrier frequencies.It is therefore an important challenge to develop high-performance equalization schemes for doubly selective channels.
Doubly selective channels can rise both in single-carrier systems and in Orthogonal Frequency Division Multiplexing (OFDM) systems.In single-carrier systems, the doubly selective channel is modeled as a time-varying filter and introduces time-varying Inter Symbol Interference (ISI).In OFDM systems, the time selectivity of the channel destroys the orthogonality between subcarriers and introduces Inter Carrier Interference (ICI) while the frequency selectivity of the channel causes the ICI to be frequency varying.In this paper, we concentrate on single-carrier systems only.
The problem of equalization for doubly selective channels has been extensively researched.Several methods for only training-based-equalization were proposed in [1,2].Semi-blind equalization methods, that can benefit from both the training and data symbols, were proposed based on linear processing [3,4] and Decision-Feedback Equalization (DFE) [5].However, the performance of these equalization methods may not be satisfactory, especially when only one receiving antenna is present [6].Moreover, the constant advance in processing power calls for more sophisticated equalization schemes that can increase network capacity.
Maximum-likelihood detection based on the Viterbi algorithm is a widely known technique for slowly fading channels.For higher Doppler rates, this method is not satisfactory due to the inherent delay in the Viterbi detector which causes the channel estimator part not to track the channel sufficiently fast.A partial remedy is offered by the Per Survivor Processing (PSP) approach, proposed originally in [7] and justified theoretically from the Expectation-Maximization (EM) principle in [8].Using the PSP approach, the channel estimation is updated along each survivor path each symbol period using Least Mean Square (LMS) or Recursive Least Squares (RLS) [9].However, for high Doppler rates, the performance is limited as these algorithms are not able to track fast fading channels [10].Improved performance can be gained by using Kalman filtering [9,[11][12][13] but this approach requires the knowledge of the channel statistics which is normally not known a priori and its estimation will likely not be able to track fast fading.Low-complexity alternatives to the PSP were also proposed [11,14].
Fast time-varying channel estimation might be achieved using the basis expansion (BE) model [15].In this method, the channel's time behavior is modeled as a linear combination of basis functions.Basis functions can be polynomials [15], oversampled complex exponentials [2,16], discrete prolate spheroidal sequences [17], and Karhunen-Loeve decomposition of the fading correlation matrix [10].Several receiver structures were proposed based on the combination of the BE with Viterbi algorithm variants like PSP [10], M-algorithm [14], and minimum survivor sequence [6,18].One common drawback of these methods is that the channel estimation part uses only hard decisions and does not weight the probability of different hypothesis for the symbol sequences.Moreover, all of these methods provide only hard decisions outputs which make them unsuitable for coded systems.
In order to enable the channel estimation part to benefit from soft decisions, several MAP-based algorithms combined with recursive, RLS-based, channel estimation were proposed [19][20][21][22].The combination of MAP decoding and maximum likelihood channel estimation can be justified using the EM principle.This leads to an iterative detection and channel estimation algorithm based on the Baum-Welch (BW) algorithm, proposed in [23] and modified for reduced complexity in [24] for non-time-varying channels (See [25] and references therein for more non-time-varying semiblind equalization methods).Adaptation for doubly selective channels is found in [26] based on incorporation of LMS and RLS in the algorithm.
Iterative MAP detection combined with polynomial BE was proposed in [27].Unfortunately, this method cannot be directly extended to higher order BE models required in high mobility environments because the choice of polynomial expansion creates numerical difficulties for higher BE models.Furthermore, the equalization and channel estimation in [27] are done in a two-step ad hoc approach which is not a true EM (see Appendix A) and exhibits degraded performance in our simulations.
Finally, Turbo equalization schemes, encompassing iterative detection and decoding, were proposed based on RLS/LMS channel estimation [22,26] and BE channel estimation [28].The latter method employs a low-complexity approximation to the MAP algorithm for the detection part.It requires, however, that the channel statistics is fully known a priori.
In this paper, we present a novel method for semiblind ML-based joint channel estimation and equalization for doubly selective channels.The method is based on an adaptation of the EM-based algorithm for doubly selective channels by incorporating a BE model of the channel in the EM iterations.Using the BE method, we can simultaneously use long blocks thereby enhancing the performance in noisy environments without compromising the ability to track the channel because of the usage of sufficiently high-order BE to model the channel time variations.The proposed method is shown to have superior performance over previously proposed methods with the same block size and number of pilots in the block.Alternatively, it requires a lower number of pilots to achieve the same performance thereby enabling more bandwidth for the information.In addition, it is shown to have good performance for relatively small blocks, which is important if low latency in the communication system is required.The proposed algorithm outputs are the log-likelihood ratios (LLRs) of the transmitted bits, making it ideally suited for coded systems and also suitable as a building block for Turbo equalization algorithm that iterates between detection and decoding stages to improve the performance further.We treat the case of uncorrelated channels paths which is the worst case in terms of number of required BE functions.In Appendix B, we discuss the generalization to correlated paths.
Another contribution of the paper is the determination of a pilot positioning scheme that improves the equalizer's performance.In the context of our proposed algorithm, the main purpose of the pilots is the enablement of sufficient quality initialization of the EM iterations so that the probability of convergence to a local maximum is minimized.To that end, we propose an initialization scheme based on a small number of pilots and find the optimal pilot positioning such that the initial channel parameters guess is as close as possible to the channel parameters obtained assuming perfect knowledge of transmitted symbols (this is the channel estimation expected at the end of the BW iterations).It is shown that the pilot positioning depends on the channel's Doppler.For high Doppler rates, our results indicate that spreading the pilot symbols evenly throughout the block leads to the best initial channel guess.This result is surprising as it is different from previous results where the optimal positioning scheme was found to be spreading of groups of pilots whose length depended on the channels delay spread [29,30].These previous results, however, were obtained using different criteria and channel model.More importantly, the analysis in these papers was restricted and did not consider pilot groups shorter than the channel's delay spread as done in this paper.Therefore, these previous results do no contradict with our new result.
Pilot positioning was discussed in [31][32][33] and conditions for MMSE optimality of both the pilot sequence and positioning were derived.The resulting sequences and positioning are, however, less attractive for practical implementations.This is because most of them require that the pilots and data overlap in time which complicated the receiver structure.The only optimal scheme proposed in [32,33] with nonoverlapping pilots and data requires a pilot pattern that results in very high peak to average transmission which is not desirable in practical communications systems.In contrast, we optimize the pilot positioning given a predefined pilot sequence (in this paper we use, as an example, Barker sequence).This allows us to derive optimal pilot positioning for a given pilot sequence that meets some other constrains (e.g., constant envelope signals, low peakto-average ratio, etc.).
The rest of the paper is organized as follows.In Section 2, we present the system model and introduce the BE model.In Section 3 we present our proposed method for semiblind joint channel estimation and equalization for doubly selective channels.Section 4 discussed the BE functions set selection.Our results regarding optimal pilot placement are presented in Section 5. Section 6 presents our simulation results and conclusions are drawn in Section 7. Partial results of this work were introduced in a conference paper [34].

Problem Formulation
2.1.System Model.The transmitted symbols vector x = [x 0 , . . ., x N−1 ] T is an i.i.d sequence with uniform distribution over an arbitrary constellation of size M.The sequence is transmitted over an unknown multipath channel modeled as a time-varying finite impulse response (FIR) filter with coefficients vector at time sample n, The received sample at time n is where y = [y 0 , . . ., y N−1 ] T is the received vector (observation vector) and x n = [x n , . . ., x n−L+1 ] T represents a branch (transition) on the trellis formed by the channel's memory [23].There are M L possible branches at each time sample n.Each possible branch is denoted by the row vector s k,n , where 0 ≤ k < M L and 0 ≤ n < N. Finally, w = [w 0 , . . ., w N−1 ] T is an Additive White Gaussian Noise (AWGN) sequence with zero mean and an unknown variance σ 2 .The time selectivity of the channel is typically characterized by the normalized Doppler frequency defined as where f c is the communication system carrier frequency, v is the user's velocity, c is the speed of light, and T s is the time of one symbol.The sequence h (i) = [h i,0 , . . ., h i,N−1 ] T , which represents the time variations of the ith channel's path, is modeled as a wide-sense stationary stochastic process with autocorrelation function [35] where J 0 is the zero-order Bessel function and Δn is the time difference in sample units.Furthermore, α i is the average power of the ith channel path and the power profile of the channel is α = [α 0 , . . ., α L−1 ] T .In addition, we make the following standard assumptions.
(A1) Information symbols, channel realization, and noise samples are statistically independent.

Basis Expansion Model.
Using the BE approach, we model the time variation of each channel's path with a linear combination of several basis functions, that is, the value of the ith path at time n is where b n (q) is the nth element from the qth basis and g i,q is the combination coefficient of the ith path and the qth basis function.The advantage of this description is that the complete time and frequency behavior of the channel is described using a relatively small set of LQ coefficients vector g = [g 0,0 , . . ., g 0,Q−1 , g 1,0 , . . ., g L−1,Q−1 ] T .We further define a BE matrix as and ] is a row vector of the function values at time n.Equation ( 4) in matrix form is then where

The Baum-Welch Algorithm for Equalization of Doubly Selective Channels
In this section, we present our new algorithm for semiblind maximum-likelihood joint channel estimation and equalization for doubly selective channels.We treat channels with uncorrelated paths as this is the worst case in term of number of required basis functions in the BE description.
In Appendix B, we extend the algorithm for channels with correlated paths.

Algorithm for Blind Equalization.
If we define the ML estimation of the channel parameters as θ = [ g T , σ 2 ] T , we would like to find θ such that is maximized.The sum is over all possible transmitted symbols vectors, or equivalently over all possible transition sequences in the trellis S. Direct maximization of p(y | θ) is an intractable problem.We can, however, maximize this expression iteratively using the EM algorithm [23].In each iteration we compute, in the E step, the expectation of the log-likelihood of the complete data conditioned on EURASIP Journal on Advances in Signal Processing the observation and our current estimate of θ.At the lth iteration, this value can be shown to be [23] Q We may express h n as where I L is an L × L identity matrix and the sign "⊗" represents a Kronecker product.In the M step, we find new θ such that this expression is maximized, that is where l is the iteration index.We may now use the new time-varying expression for the channel to get the doubly selective version of the algorithm in [23].Plugging ( 9) in ( 8), repeating the derivation in [23], and utilizing the Kronecker product properties, the resulting update equations are where we have used the identity xy ⊗ zw = (x ⊗ z)(y ⊗ w) and The values of p(s k,n | y, θ (l) ) are efficiently computed using the Bahl-Cocke-Jelinek-Raviv (BCJR) algorithm [36].For non-time-varying channels, we have b n = 1 and (11) reduces to equation (11) in [23].Our algorithm may also be extended to channels with correlated paths.In this case, the number of BE parameters used for describing all channel's paths may be reduced.In that sense, uncorrelated channels paths may be considered as the worst case.The correlated case is discussed in Appendix B.

Adaptation to the Semiblind
Case.Adaptation of the above algorithm to the semi-blind case where we have some known pilot symbols is straightforward.The only required change is in the computation of p(s k,n | y, θ (l) ) using the BCJR algorithm.We modify the branch metrics so that all transitions that are not consistent with known pilots are assigned zero probability.This ensures that transition probabilities are calculated with the a priori information about the pilots.

Initialization of the Algorithm.
Optimization of the EM objective function ( 8) is a nonlinear process that may converge to a local maximum.It is therefore important to calculate good initial guess for the channel parameters so that the probability of convergence to a local maximum is minimized.We suggest using the available pilot symbols for finding initial channel parameters using the following method.First, we run the BCJR algorithm where the branch metrics are initialized without any initial channel guess by assigning zero a priori probability to all transitions that are not consistent with the known pilots and equal (nonzero) a priori probability to all transitions that are consistent with the pilots.This initialization of the branch metrics represents our best a priori knowledge about the transitions probability in the trellis.From the BCJR algorithm we find p(s k,n ) and then use them in (11) to obtain an initial guess for the channel BE parameters.More details on this initialization method can be found in Appendix C.An important feature of this initialization scheme is that all observations that have some content of pilots in them are taken into account including those with mixed pilot and data contributions.This is in contrast to most other pilot based estimations that take into account observations based on pilots only [30].This fact turns out to be significant when we discuss how to position the pilot symbols in the block in Section 5.
It should be noted that the above algorithm requires initial synchronization stage to ensure that all major channels taps fall within the searched multipath window.This synchronization stage can be done at a much lower rate than channel estimation update, as the channel tap positions typically drift at a much slower rate compared to the fading rate, and therefore its complexity is negligible.The synchronizations stage is outside the scope of this paper and we assume perfect synchronization throughout the paper.

Computational Complexity.
The computational complexity of the updating equations ( 11) and ( 12) is analyzed in Table 1, where we have broken the calculation to several stages and counted the number of complex Multiply-And-Add operations (MAC) for each stage.For comparison, the equivalent complexity of [23] can be obtained from the same table by eliminating the stages for calculating T 2 , T 4 , and h n .For (11), it can be seen that for the typical case of M L > Q 2 /2, the stages of calculating T 1,n and T 3,n are more computationally complex than the stages of calculating T 2 and T 4 , respectively.For (12), it can be seen that the second stage of calculating σ 2 is more complex than the first stage of calculating h n .Finally, note that all the computationally complex stages of T 1,n , T 3,n , and σ 2 do not depend on Q and therefore their complexity is the same as in [23].We may therefore conclude that the proposed algorithm extends the Baum-Welch algorithm [23] for doubly selective channels with only minor increase in complexity.

Selection of the Basis Functions
Although many basis functions are possible, previous papers concentrated mostly on three types of basis function sets.The first one is the complex exponentials functions set [2,37].The value of the qth basis function of this set at time n is These functions are periodic with period N bem .In order to avoid modeling errors at the block edges, we therefore set N bem = 2N.The second type of basis functions is [15] b n q = n q .( The functions in (14) model the channel time behavior as polynomial in time.This choice of basis functions may be regarded as a generalization of the channel description in [27] where it was suggested to use first-and second-order polynomials to model the channel time variations.The best basis functions are the ones that minimizes the mean square error of the fading process description given a finite set of Q basis functions.That is, where the vector h represents the channel time variation.The solution for this problem is readily available by usage of the Karhunen-Loeve Transform (KLT) [10] and the basis functions are the eigenvectors of the autocorrelation matrix of the Rayleigh fading process.The element n1, n2 in the autocorrelation matrix is Out of all eigenvectors, the Q vectors that correspond to the largest eigenvalues are selected as the basis set.The target function in (15) is suitable for flat fading channel.
For frequency-selective channel, the mean square error will be simply the sum of the mean square errors of the individual paths and therefore the same solution is optimal for multipath channels.We note that a similar argument is given in [38].
An obvious alternative to the equalization approach we propose in this paper is to divide the data block into small subblocks such that the channel can be considered approximately constant within a subblock period and then equalize each subblock separately using the Baum-Welch algorithm for non-time-varying channels [23].Interestingly, this subblock scheme can be considered as an instance of the BE approach if we choose Q basis functions for Q subblocks where the qth basis function is equal to one in the symbols time that correspond to the qth subblock and zero elsewhere, that is , where 1 x is a vector on x ones.To justify our approach, we would like to compare it to this subblock approach.
The efficiency of a given set of basis functions may be evaluated by calculating the mean square error of the fading process representation using this set of functions: where E and Tr are expectation and matrix trace operators, respectively.Figure 1 plots the required number of basis functions (rank of B) so that the mean square error in (17) is lower than 1% error.As expected, using the eigenvectors as basis functions leads to the lowest number of functions.
The polynomial basis set is shown to be quite close to the optimal eigenvectors solution for low normalized Doppler while for high Doppler rates, it is more beneficial to use the complex exponentials basis.The sub-block-based basis functions performance is much worse.This is not surprising as these basis functions do not utilize the correlation between subblocks and force a noncontinuous description of the channel in contrast to the channel's typical behavior.The results shown in Figure 1 confirm that this choice of basis functions is not suitable for Rayleigh fading and provides an explanation to the degraded performance of the subblock method shown in the simulation results section.probability of convergence to a local maximum is minimized.First, we reformulate the initialization scheme in Section 3.3 as an equivalent Least Squares (LS) problem.Consider first the case where all transmitted symbols are known.In this case, the channel parameters g can be found with an LS solution to the problem

Placement of Pilot Symbols
where A represents the known transmitted symbols and the BE model-based time variations.More specifically, where X = [X 0 , . . ., X L−1 ] and X n is an N ×N diagonal matrix such that X n = diag(x −n , x −n+1 , . . ., x N−1−n ) and negative indexes represent data symbols from the previous block that are affecting the observations of the current block due to ISI (if no interblock interference is assumed, these can be replaced by zeros).In addition, B L = I L×L ⊗ B and I L×L is an L × L identity matrix.The LS solution to ( 18) is Now consider the case where only part of the transmitted symbols are known (pilots) and replace the unknown symbols in X with zeros.The solution to this "sparse LS" problem is where and X p is defined similarly to X with nonpilot symbols set to zero.Finally, B p is received by setting to zero all elements in the rows corresponding to nonpilot symbols in the matrix B L .In Appendix C, we show that the initialization method in Section 3.3 is equivalent to (21).The initialization method is thus equivalent to finding the best BE model parameters vector that fits, in the LS sense, the transmitted pilot sequence (Note that the noise term in this model is not white (since the data is treated as part of the noise).Therefore, a better initialization would be to use weighted least squares method.
To do that, however, the noise level and average channel profile need to be known or estimated).Our goal is to position the pilots such that the initial channel guess, based on these pilots, will be optimal according to some criterion.Two reasonable criteria for pilot positioning are where p is a vector of the pilot positions in the block.The maximum function in (23) and expectation in (24) are taken with respect to the data symbols, noise, and channel realizations.Using these criteria, it might be possible to optimize both the pilots positions and the pilot patterns.We, however, select known pilot patterns (e.g., Barker sequences) so that we keep constant envelope signals and optimize the positioning for this given pilot pattern.The usage of these two criteria is detailed in the next sections.Interestingly, both criteria lead to the same positioning scheme for high Doppler rates.
EURASIP Journal on Advances in Signal Processing 7

Worst Case Analysis for Flat Fading Channels.
In this section, we find the best positioning scheme by using (23).First, notice that the criterion may be decomposed to two terms because where the second equality is justified because, by construction of g, the term y − Ag is orthogonal to the span of the matrix A to which Ag − Ag p belongs.Note that only the second term is dependent on the pilot positions.Obviously, the best pilot positioning is dependent on the channel and noise realizations.Our goal is to obtain positioning scheme suitable for all channels, data, and noise realizations by optimizing the positioning scheme with respect to the worst case realizations.Using ( 19) and ( 22), the second term in ( 25) may be bounded by where and σ 2 max is the largest eigenvalue of the matrix A H p A p .For flat fading channels (and any PSK constellation) X H X = I and, therefore, where eig[D] is the vector of eigenvalues of the matrix D.
The second equality follows from the fact that for flat fading channels Matrices D and X −1 DX have therefore the same eigenvalues).It follows that minimization of the worst case MSE is achieved by finding a pilot positions vector p such that The matrix B H p B p is a deterministic function of the BE functions, block size, pilot positioning, and the pilot pattern (sequence).It is therefore possible to find the best positioning scheme for the desired block size, BE model, and pilot sequence with a computer search.For simplicity, we limit the search for patterns in which the pilots are grouped in groups of length L and these groups are spread throughout the block as evenly as possible.This means that the pilot positioning we find with this limited search is only optimal amongst all positioning with evenly spaced pilot clusters.However, all previous works on pilot positioning arrived at positioning schemes that are consistent with this structure.It turns out that the best positioning scheme is obtained with     L = 1 for all tested block sizes.It is interesting to note that this result is identical to the result in [29] which was obtained using different channel model and criterion.

Mean Case Analysis for Frequency Selective and Frequency
Flat Channels.In this section, we optimize (24).We begin with the approximation where and x p (m) is defined as     Note that an accurate expression (with no approximation) may be obtained by replacing x * p (n − k)x p (n − j) with x * n−k x n− j .When either x n−k or x n− j is an information symbol (not a pilot), this multiplication result is a random variable, uniformly distributed over a finite set of values with zero average.As a result, for long enough blocks, the contributions from the information symbols to the sum in (30) cancel out and this approximation is fairly accurate.Our criterion may be therefore approximated with The analysis that follows should be considered valid only for large enough block sizes where ( 29) is accurate.The expectation of the approximated metric is The autocorrelation matrix k yy H X j .Using the standard assumption that the channel's paths are statistically independent (assumption A2), we may express the autocorrelation matrix R as a linear combination of the contributions of the channel paths, that is, Using assumptions A1-A2 and (3), the entry n1, n2 in the submatrix k, j (or equivalently, the element kN + n1, jN where The best pilot positioning scheme is therefore This expression is deterministic and depends only on the BE functions, block size, noise variance, channel order (L), Doppler rate, and pilot sequence.It is therefore possible to find the best pilot positioning for a particular set of parameters by evaluating (37) for various p.As we did in the previous section, we limit the positioning patterns for patterns in which the pilots are grouped in groups of length L, and these groups are spread evenly throughout the block.This positioning strategy coincides with the pilot positioning in [39] for L = 1, with the pilot positioning in [30] for L = L and with the pilot positioning in [29] for L = 2L + 1.In addition, every group of pilots is a Barker sequence of length L. Barker sequences are known to enable good channel estimation because of their autocorrelation properties.Define the positioning metric as the value to be minimized in (37).A typical behavior of this positioning metric is shown in Figure 2 (based on (37) and in agreement with simulation results).
The optimal positioning strategy is shown to be dependent on the Doppler rate and the number of pilots in the block.As can be seen from Figure 2, for low Doppler rates it is better to use group of pilots as also indicated by [29,30] (although the difference is not very significant, at least for short delay spreads).For high Doppler rates and a small number of pilots, however, it turns out that using L = 1 leads to much better results.This is because there is a tradeoff between accurate estimation of the multipath at specific points in time (that is better achieved by grouping the pilots) and tracking the channel time variations (that is better achieved by spreading the pilots throughout the block).Our results indicate that for high velocities using L = 1 leads to a lower metric value as this means better ability to track time variations.Note that this result is obtained for severe ISI channel with three equal energy paths (and similar result was obtained for channel with 5 equal energy paths).We have also simulated channels with less severe ISI (that is, decaying power profiles), and the advantage of using L = 1 was even larger, as could be expected.The switching point (Doppler rate beyond which it is advantageous to use L = 1) is dependent mainly on the percentage of pilots in the block.For larger number of pilots, the switching point will occur at higher Doppler rate.The reason is that for large number of pilots there will be sufficient number of groups in the block to allow tracking of path time variations even when the group size is kept 2L + 1, so both multipath profile and time variations could be estimated accurately.We, however, are interested in the smallest number of pilots that enables good performance, and in these conditions, L = 1 is advantageous even for moderate Doppler rates (see Figure 2).This conclusion is somewhat surprising as it is different from previous conclusions in [29,30].However, these previous works used different channel models and performance criteria.Moreover, both works considered only pilot groups equal to 2L + 1 [29] or L [30] or longer, to facilitate their analysis.

Simulation Results
6.1.Performance of the Proposed Equalization Scheme.Next, we present simulation results for our proposed equalization scheme.We use a sequence of 2 17 QPSK symbols that is sent through a doubly selective channel as described in Section 2.1.The normalized Doppler frequency is f nd = 0.002, and coherence time, defined as the time over which the channels response to a sinusoid, has a correlation greater than 0.5 is 9/(16π f nd ) = 96 symbols.A modified Jakes fading model is used to model the time variations of each of the channel paths [40].The pilots are positioned according to the optimal scheme found in the previous section (L = 1).The number of basis functions for all simulated BE sets is where N bem = 2N.This number was tested numerically to enable good accuracy description of the channel with the BE complex exponents and polynomial functions (below 1% error).This is also the number of basis functions used in [30].For the selection of eigenvectors as the functions set, we could have decreased this number slightly.We present simulation results for the following equalization algorithms.
(i) Maximum Likelihood equalization using perfect channel knowledge.
(ii) Maximum likelihood equalization with channel estimation based only on the pilots.This is identical to the first iteration of the proposed algorithm.
(iii) Time-varying BW algorithm with BE based on complex exponential functions (13) (BW-BE-exp).
(iv) Time-varying BW algorithm with BE based on complex exponential functions (13) and initial channel guess identical to the true channel (BW-BE-exp & perfect init.).The difference between the error curve of this simulation and the previous one will indicate if we have an issue of convergence to a local maximum.
(v) Time-varying BW algorithm with BE based on polynomial functions (14) (BW-BE-poly).This might be considered a significant improvement of [27].
(vi) Time-varying BW algorithm with BE based on optimal basis functions (BW-BE-eig).
(vii) Non-time-varying BW algorithm based on dividing the data blocks into shorter blocks in which channel is assumed to be constant (BW-SB).This is essentially the method of [23].
(viii) The BW-RLS method in [26] (called APP-SDD-RLS in [26]).This method was initialized using the same initialization scheme we used for the BW-BE methods.After the parameters of the BE are found, the actual channel responsed estimate is computed for every time instance.Finally, the BCJR algorithm uses this estimate to calculate the transitions probabilities which are the starting point for the BW-RLS in [26].
(x) Iterative Viterbi-based equalization with BE-based channel estimation.Reduced complexity variants of this algorithm appeared in [14,18].
Simulation results for various signal-to-noise ratios (SNR) are presented in Figure 3 for block size of 256 symbols, 20 pilots (about 8% pilots), and multipath channel with three symbol-spaced paths with power profile [0, −3, −3] dB.The proposed BE-based EM algorithm performance is very close to the performance of a nonblind ML equalizer with perfect channel knowledge for both complex exponentials and polynomial BE functions.The BW-RLS method and the BW-SB have similar performance.This is not surprising as the BW-RLS method can be considered a non-timevarying BW that operates on a (exponentially weighted) sliding block whose length is proportional to the forgetting factor used.The BW-RLS can be therefore considered as a generalization of the BW-SB.The PSP-RLS method exhibits poor performance in these scenarios as the RLS component is not able to track the fast fading.The performance of Vit-BE is closer to our proposed methods although still significantly inferior.This is because the channel estimation part uses only hard decisions and does not weight the probability of different hypothesis for the transmitted sequence.The training based ML equalizer has the worst performance as the number of pilots is not sufficient for reliable channel estimation.It is therefore evident that the proposed semiblind approach can reduce the number of required pilots significantly thereby increasing the available bandwidth for information.
Similar results were obtained for different channel profiles and delay spreads.For example, Figures 4 and 5 present the performance of the proposed algorithm for two equal power taps and four equal power taps, respectively.It is clearly shown that the proposed algorithm outperforms all previously proposed schemes that were simulated.
Results for block size 256 and SNR = 12 dB are shown in Figure 6 for various numbers of pilot symbols (the zero pilot percentage point in the figure corresponds to totally blind case and uses random initial guess).It is shown that in these conditions, our proposed method requires about 12% pilots to achieve close to optimal performance while Vit-BE requires up to 20% pilots to reach comparable performance and the other methods performance is even worse.Trainingbased equalization does not produce reliable results because the fast fading channel necessitates an excessive number of pilots for reliable estimation and tracking of the channel throughout the block in these conditions.
Results for various block sizes are shown in Figure 7.The performance of the proposed method improves with increasing the block size even though this means more functions are needed for modeling the time behavior of the channel.The reason is that the larger block size enables better filtering out of the noise as well as benefiting from more pilots in the block.All simulated methods reach a saturation point beyond which it is not beneficial to increase the block size further.The proposed method is shown to reach this point for smaller blocks thereby enabling latency reduction.From other simulation results (not shown here as they resemble the results presented so far), it is concluded that the proposed methods of BW-BE are superior over all simulated previously proposed methods in a large range of block sizes, pilot percentages, and channel profiles.
Finally, although a detailed convergence analysis is beyond the scope of this paper, a histogram of the number of iterations required to achieve convergence is plotted in Figure 8.This plot is for SNR = 12 dB and under the same conditions as in Figure 3. Convergence is determined based on the Q function value, defined in (8).Convergence is assumed when the ratio of two Q function values of two consecutive iterations is sufficiently close to one, that is, less than 1 + 10 −10 .The maximum number of iterations allowed in this plot is 30.It is shown that about 75% of the blocks converge with 8 iterations or less and 90% of the blocks require 15 iterations or less for convergence.These results can be used to evaluate the computational complexity of the proposed algorithm.

Conclusion
A novel method for maximum likelihood semi-blind joint channel estimation and equalization for doubly selective channels is proposed.This method is based on expectation maximization (EM) principle combined with a BE model of the doubly selective channel.It is shown that the main drawback of the proposed approach is its possible convergence to local maxima.To alleviate this problem, we propose an initialization scheme to the equalization based on a small number of pilot symbols.We discuss the selection of basis functions set and find the optimal basis set.We further derive a pilot positioning scheme targeted to reduce the probability of convergence to a local maximum.The resulting algorithm is shown to be significantly superior over previously proposed equalization schemes and to perform in many cases close to an ML equalizer with perfect channel knowledge.

A. Comparison of the proposed method to [27]
In this appendix, we highlight the differences between the algorithm proposed in [27] and the algorithm presented in this paper.For simplicity, we discuss only the non-timevarying case as the generalization to time-varying case is straightforward.
The channel estimation is updated in the BW algorithm at each iteration with: In [27], the channel is updated in two stages.First, the expectations of the Hidden Markov Model states are computed with where m k is the expected received value when in state k.
In the second stage, the channel coefficients are updated by Comparing (A.5) to (A.1), we see that the two methods are not equivalent.Moreover, the method in [27] converges to the true BW only when N−1 n=0 p(s k,n | y, θ (l) ) is constant for all values of k.This is likely to happen for long blocks.However, for short block it is not necessarily true and in this case the approach presented in [27] leads to a suboptimal solution.

B. Extension to correlated channel's paths
In this appendix, we extend the proposed equalization algorithm for the case of correlated channel's paths.First, let us assume for simplicity we have only two paths in the channel and rewrite (6)  where now h = [h (0)T , h (1)T , . . ., h (L−1)T ] T and h (i) is a a vector of size N that represents the ith channel tap amplitude time variations.The size of the matrix B now is LN × Q.
The main difference in the correlated case compared to the uncorrelated case is that the matrix B is no longer block diagonal and the off diagonal elements represent the correlation between the taps.In this description, there is a common coefficients set used for all channel taps and the difference between taps is due to different basis functions for each tap (the basis functions for the ith tap are the rows iN to (i+1)N −1 in the matrix B).The criterion for basis selection, similar to (15), is Obviously, Q needs to satisfy Q ≤ Q ≤ QL and will vary according to how much the taps are actually correlated.
In the uncorrelated case, we will have Q = QL and in fully correlated case Q = Q.With uncorrelated taps, the autocorrelation matrix E[hh H ] is block diagonal and there are L blocks along the diagonal, each is equal to R corr as defined in (5).The eigenvectors matrix in this case will also form a block diagonal matrix with L identical blocks whose first Q rows are the basis functions of the individual taps.In this case, the basis functions for each tap are identical and we return to the description of the channel response given in (9).In the correlated case, we have

C. Proof of the proposed BW initialization method
In this appendix, we prove the equivalence between the initialization method of the BW based on symbols priors, discussed in Section 3.3 and the initialization method shown in (21).First, define the group n = {n, n−1, . . ., n−L+1}, that is, the time instances on which the transition s k,n depends.In addition, define the group G as the group of all transitions in the trellis that are consistent with the pilots.Our branch metric initialization is The probability of a transition in the trellis is the sum of probabilities of all paths going trough this transition.The probability of every path that is consistent with the pilots is

Figure 2 :
Figure 2: Pilot positioning metric (the value to be minimized in (37)) for various L, block size = 512, 5% pilots in the block, channel order L = 3, equal average energy paths, number of basis functions Q = 2 N bem f nd + 1.
-exp and perfect init.
-exp and perfect init.
BW-BE-exp and perfect init.
BW-BE-exp and perfect init.

12 EURASIP
Journal on Advances in Signal Processing fitting them to the values of m = [m 0 , . . ., m M L −1 ] T in a least squares sense, that is, h = arg min h |Sh − m| 2 , (A.3) where S is an M L × L matrix whose rows are the states of the trellis, that is, S = [s T 0 , . . ., s T M L −1 ] T .The well-known solution is S H S h = S H m. (A.4) Using the definition of S and plugging (A.2) into (A.4),we obtain k,n | y, θ (l) y n N−1 n=0 p s k,n | y, θ (l) .(A.5)

2 .
b i is the ith row in the matrix B. The matrix B n is then built from the L rows in B that correspond to time instant n.This description can then be used in the derivation of the EM procedure and the resulting updating equations areg (l+1) ,n | y, θ (l) B H n s H k,n s k,n B n ,n | y, θ (l) y n − s k,n B n g C N and p s k,n = ⎧ ⎨ ⎩ M N−|p n| C N , s k,n ∈ G,
Figure 1: Required number of basis functions for mean square error less than 1%.Block size = 256.
[29,30,39]pread evenly throughout the data block is a better strategy[29,30,39].Proper pilot placement for EM based algorithms is particularly important because of the highly nonlinear nature of the EM objective function in doubly selective channels, which results in many local maxima.The purpose of the pilots is, therefore, to enable sufficient quality channel parameters vector initialization so that the