Group-Sparse Channel Estimation using Bayesian Matching Pursuit for OFDM Systems

Abstract

We apply the Bayesian matching pursuit (BMP) algorithm to the estimation of time-frequency selective channels in orthogonal frequency division multiplexing (OFDM) systems. By exploiting prior statistics and sparse characteristics of propagation channels, the Bayesian method provides a more accurate and efficient detection of the channel status information (CSI) than do conventional sparse channel estimation methods that are based on compressive sensing (CS) technologies. Using a reasonable approximation of the system model and a skillfully designed pilot arrangement, the proposed estimation scheme is able to address the Doppler-induced inter-carrier interference (ICI) with a relatively low complexity. Moreover, to further reduce the computational cost of the channel estimation, we make some modifications to the BMP algorithm. The modified algorithm can make good use of the group-sparse structure of doubly selective channels and thus reconstruct the CSI more efficiently than does the original BMP algorithm, which treats the sparse signals in the conventional manner and ignores the specific structure of their sparsity patterns. Numerical results demonstrate that the proposed Bayesian estimation has a good performance over rapidly time-varying channels.

1. Introduction

Channel estimation for high mobility wireless communication scenarios is a challenging task in orthogonal frequency division multiplexing (OFDM) systems because the severe Doppler frequency shift destroys the orthogonality of subcarriers, which causes strong inter-carrier interference (ICI) to pilot symbols. This interference results in a sharp decline of the estimation performance. To address the ICI, traditional estimation methods [1] such as the least square (LS) method require a large number of pilot symbols, thereby leading to a reduced bandwidth efficiency.

Recently, the compressive sensing (CS [2, 3])-based channel estimation [4-6] has gained a fast-growing interest. Compared with traditional pilot-aided methods, the CS-based channel estimation is considered as a more promising technology because it is able to take advantage of the inherent sparsity of the propagation channel. This superiority enables the latter to detect the channel status information (CSI) with considerably reduced pilot symbols in high mobility communication environments. Moreover, in the recent work [7, 8], experimental results demonstrated that CS-based methods are clearly superior to traditional estimation methods over some real-world multipath channels (e.g., the underwater acoustic channel) in OFDM systems.

As an extension of CS, the Bayesian CS (BCS) [9-11] technologies such as the Bayesian matching pursuit (BMP) [12-14] are based on prior statistics of the sparse signal. Using the additional priori statistical knowledge, the Bayesian recovery method outperforms the conventional sparse recovery methods in many applications. In the field of channel estimation, some statistical information of the channel can also be obtained in some scenarios, providing the potential for a more precise estimation of the channel or a further reduction of pilot symbols. Based on the available channel statistics, the recent work presented in [6] developed a statistical basis optimization procedure, which enables channel coefficients to be represented by a sparser basis expansion model. Simulation results in this literature showed that the statistically optimized basis can achieve an improvement of the estimation performance for doubly selective channels. However, until now there has been little CS-based work that considers the direct application of BCS approaches to the channel estimation. Therefore, we believe that the statistical information of the channel is not effectively exploited in the existing sparse channel estimation methods.

Compared with the related work mentioned above, the uniqueness of the approach presented in this paper can be summarized as follows. First, using the channel basis expansion function in the delay-Doppler domain, we depict the group-sparsity (i.e., block sparsity) [15, 16] and statistical properties of doubly selective channels. Based on these priors, we represent the CSI by a statistical group-sparsity model, which enables us to formulate the channel estimation as a problem of reconstructing the statistical group-sparsity signal. Then, we propose a BMP-based method to reconstruct the CSI for rapidly time-varying OFDM systems. In this estimation scheme, we adopt an approximate model of the fast fading transmission system. This approximation allows us to handle the Doppler-induced ICI with a relatively low computational complexity. Moreover, because the BMP algorithm treats the signal as being sparse in the conventional sense and ignores the specific structure of their sparsity patterns, we make some small changes to BMP, which enable the algorithm to make good use of the group-sparse structure of the channel. Simulation results demonstrate that these modifications lead to a further reduction of the computational cost. Finally, we also analyze the influence of different pilot patterns on the estimation performance and then present an approximately optimal pilot placement for the proposed estimation scheme.

The remainder of the paper is organized as follows. We introduce some fundamentals about the BMP algorithm and the OFDM system in the next section. The sparse and statistical properties of the doubly selective channel are described in Section 3. In Section 4, we propose a low-complexity estimation method which is based on the modified BMP algorithm. The pilot design is also studied. Section 5 shows the experimental results. Conclusions are finally drawn in Section 6.

 

2. Preliminaries

2.1 BMP Method

CS and other sparse recovery methods focus on the classical linear measurement model defined as follows

where u ∈ ℂM is an unknown K-sparse vector (has only K nonzero entries, K << M), Y ∈ ℂN ( K < N < M ) is a known measurement signal, Փ ∈ ℂN×M is the measurement matrix and is an additive Gaussian noise. It has been proved in [2] that u can be reconstructed from Y by conventional sparse recovery methods [2, 3] if the measurement matrix satisfies the so-called restricted isometry property (RIP) [17].

In some scenarios, the sparse signal exhibits an additional structure in the form that the nonzero elements occur in groups, not arbitrary places. By taking advantage of the structured sparsity, the group-sparse recovery methods [15, 16] can further compress the information of sparse signals; thus yield an improved recovery performance compared with traditional CS algorithms.

As another special case of CS technologies, the BMP algorithm [12-14] and other BCS methods are based on the additional statistical information of the sparse signal. Generally, u is considered as a mixed Bernoulli-Gaussian process in these approaches. More specifically, let S be the support of u, then the elements of S are assumed to be sampled from {1,2, , M} according to the Bernoulli distribution with success probability P. Under this assumption, we obtain the probability of S :

It is observed that the sparsity of u is controlled by P ( E[K] = PM ). Defining uS to be the subvector of u, which is indexed by S , we further assume that the entries of uS are drawn from an Gaussian process with zero-mean and known variance σ2. Then, the minimum mean-square error (MMSE) estimate of u from Y is given by

and the maximum a posteriori (MAP) estimate is expressed as

where

To calculate the two estimates, the expectation and the posterior must be evaluated. We first consider the calculation of the former. Because the relationship between u and Y is linear, the linear MMSE (LMMSE) estimate can be used to evaluate the expectation [14]:

where ՓS is the submatrix that consists of the columns of Փ , indexed by S , INis an N × N identity matrix, and (·)H denotes the conjugate transpose matrix.

Next, using the Bayesian rule, we have

where the factor p(Y) can be ignored since it is common for each potential support, and p(S) is obtained from (2). Then, we focus on the prior probability . Because the entries of uS are drawn from 2 CN (0,σ2), we have that is CN (0, ΣS), where

So,

Note that if the dimension of u is large, there are quite many possible supports. In this case, computing the posterior and expectation for all these supports requires an unacceptable computational cost. To simplify the estimation, the computation should be restricted to the dominant supports, i.e., the supports which have significant posteriors. Defining the dominant support set as S* , we approximate the MMSE estimate as follows

which is called the AMMSE estimate. Note that the MAP estimate can be treated as an extreme case of AMMSE, where the size of S* is reduced to one. Therefore, we only analyze the AMMSE estimate in the rest of the paper. The search procedure of S* is presented in Section 4.2.

2.2 OFDM System

We focus on the discrete time-varying transmission model

where h[n, l] denotes the channel impulse response (L is the maximum multipath delay), x[n] is the input signal, y[n] is the output signal and z[n] is the Gaussian noise with zero-mean and variance . Over a block of N symbols, the linear system can be written in matrix notation as

where the vectors x = (x[0],...,x[N-1])T , y = (y[0],...,y[N-1])T , z = (z[0],...,z[N-1])T and represents the channel matrix:

We construct an OFDM system with N subcarriers. Let be the standard discrete Fourier transform (DFT) matrix , and X, Y be the transmitted and received data vectors, respectively. Then, we have x = QHX and Y = Qy. To remove the inter-symbol interference (ISI) from the OFDM system, we assume that the length of the cyclic-prefix (CP) is longer than L. Then, over a block of an OFDM symbol, we can express the overall system model as follows

where the noise vector is Z = Qz and the overall system matrix is G = QHQH. The system model can also be written in the following form

where the non-diagonal entries of G represent the ICI response. Inserting G = QHQH into (15), we have

We define the impulse response vector as follows

and an N × NL matrix A as follows

where n = (q-1)mod(N) , and l = (q-1-n) / N. Then, (16) can be expressed in matrix notation as

The system model is changed as a function of h, which is convenient for us to address the channel estimation problem.

 

3. Sparsity and Statistical Characteristics

We construct a doubly selective channel which consists of S propagation paths. The channel impulse response is expressed as

where δ(·) is the Dirac function, αs , τs and fs denote the magnitude, multipath delay and Doppler frequency, respectively. Assuming that the sampling period is T, we approximately consider that τs/T are integers for s = 1,...,S. In this case, the discrete impulse response is given by

We first investigate the sparse structure of the channel. To describe the sparsity in the delay-Doppler domain, we use the delay-Doppler basis expansion:

where

Now we focus on the following function

which describes the energy distribution of us[l,d] in the frequency domain. Generally, the normalized Doppler frequency is |fsNT| << 1; thus the frequency offset is constrained in a small band around the zero point [6-8, 18]. An example of Es[d] is shown in Fig. 1.

Fig. 1.An example of Es[d], N = 1024, fsNT = 0.3

As observed in Fig. 1, the energy of Es[d] is concentrated in the subspace d ∈ [-D,D] for an appropriately chosen D ( [-D,-1] is the circular shift of [N - D, N -1] ). Because δ(l-τs/T)=0 if l≠τs/T, the dominant components of us[l,d] are mainly concentrated in the group {us[τs/T,0],...,us[τs/T,D], us[τs/T,N - D],...,us[τs/T,N -1]}. For all S paths, the significant elements of u[l,d] cluster in S such groups. Ignoring the small-magnitude elements, we approximately consider that u[l,d] has only S nonzero groups, which only occur in the subset [0,L-1]×([0,D]∪[N - D, N -1]) . The number of valid channel coefficients is decreased from NL to (2D+1)L.

We define a (2D+1)L-length vector as follows:

Then, only S of the groups ( u0,...,uL-1 ) are nonzero. If , u is approximately an group-sparse vector over this group partition.

Next, we investigate the channel statistical properties. Without loss of generality, we assume that for s = 1,...,S , τs and fs are uniformly sampled from {0,T,...,(L-1)T} and [-fmax , fmax] , respectively, and is a Gaussian process. Because ul is nonzero only if l = τs/T , we can approximately consider that each group of u has the probability P = S / L to be nonzero (i.e., a Bernoulli process). We assume that S is a possible support of u, which corresponds to a certain group-sparse pattern. Then, the probability of S is given by

Moreover, the entries within each nonzero group of u are Gaussian. Therefore, u can be considered as a group-version of the mixed Bernoulli-Gaussian process.

 

4. BMP-Based Channel Estimation

In this section, we propose a BMP-based estimation method to reconstruct the CSI. In the proposed estimation, we make some changes to BMP, which enable the algorithm to address the group-sparse channel coefficients in a more efficient manner. Furthermore, a discussion about the pilot pattern is also presented.

4.1 Proposed estimation scheme

We consider the representation of (19) in the delay-Doppler domain. From (22), the relationship between h and u is obtained as

where ⊗ denotes the Kronecker product, and denotes the uniform down-sampling operation from [1,N] to [1,D+1] ∪ [N - D+1,N]. Using u to replace h, we transform (19) into the form

where the N×(2D+1)L matrix Փ = A(IL ⊗ (FQ)H). We define a (2D+1)L-length vector

and, subsequently, n = (mq -1)mod(N) and l = (mq -1 - n) / N for q = 1,...,(2D+1)L . Then, the entries of Փ are expressed as follows

It is seen that X(k-n-1)mod(N)+1 ≠ Xk if n ≠ 0 . Essentially, these symbols are the direct neighbors on both sides of Xk. Hence, the corresponding entries Փk ,q describe the interference of these neighboring symbols (i.e., the ICI response) on Yk . Assume that there are NP pilots in the N-symbol block, and the pilot tones belong to the subset P = {k(1),...,k(NP)} . Selecting the measurements at the pilot positions, we form the measurement model

or

In the pth row of ՓP , the entries Փk(p),q corresponding to (k(p) - n - 1)mod(N) +1∉P (n was defined in (30)) are unknown because we know only the pilot symbols, not data symbols. These unknown entries represent the ICI from the data symbols. To construct the measurement equation without the pre-estimation of ՓP , we should ignore these unknown components. The approximate matrix is defined as follows

Using to replace ZP , we change (32) into the form

To reconstruct u from (34), we use the modified BMP algorithm, which is described in Section 4.2. With the estimate of u, the impulse response vector h is obtained via (27).

4.2 Modified BMP Algorithm

To compute the AMMSE estimate of u, we must look for appropriate supports to form the dominant set S* in (10) at first. In [12, 13], a fast greedy approach was proposed to search the best supports of different sizes. Here, to make full use of the group-sparsity of u, we make some modifications to this greedy approach. And, we will show that the modified algorithm can significantly reduce the computational complexity.

Assume that the sparsity of u is unknown. Then, all potential support sizes must be taken into consideration. Since u is a group-sparse vector, we investigate only the group-sparse samples instead of all sparse ones. Assuming the maximum possible group-sparsity of u is Smax , we consider the group-sparsity that belongs to [1, Smax]. In other words, the support sizes of u are selected from the subset {2D+1,2(2D+1),...,Smax(2D+1)}. For each potential size, we choose the support S which has the highest posterior probability. For example, when searching the optimal support of size 2D+1 (i.e., one nonzero group), we evaluate for S1 = U0 , ..., UL-1 , where Ul is the support corresponding to ul (ul was defined in (25)). The support which maximizes is considered optimal for this support size. We assume that the optimal support is S1 = Ui. When evaluating supports of size 2(2D+1) , we choose Uj ≠ Ui, which enables S2 = Ui ∪ Uj to maximize . When the support size further increases, the greedy search continues in the same manner. Finally, there are Smax supports selected to form S* .

With the dominant set, the AMMSE estimate in (10) can be calculated using the fast recursive method proposed in [12]. In contrast to AMMSE, the MAP estimate in (4) needs less steps since it considers only the optimal support instead of all dominant supports. However, the reduced complexity probably leads to performance degradation. Detailed steps of the modified algorithm are presented in Table 1.

Table 1.Modified BMP algorithm

Using the modified BMP algorithm, in each iteration of the support search procedure, we only need to calculate the posterior around L times. Compared with the modified algorithm, the original BMP algorithm requires a far larger number of calculations (about (2D+1)L times) because the latter cannot use the group-sparse structure of u. Furthermore, the iteration number of the modified BMP algorithm is also decreased from (2D+1)Smax to Smax. Hence, the total number of operations is reduced by (2D+1)2 times. Obviously, the modified BMP algorithm has a far lower computational complexity.

To improve the accuracy of the AMMSE estimate, we should expand the dominant set S* to involve more potential supports. This can be done by repeating the support search step several times. In each repetition, the supports which already belong to S* should not be chosen again. If we repeat the search procedure R times, the size of S* is increased from Smax to RSmax.

In Table 1, Smax, P, σ2 and can be estimated from some prior channel information or set as predefined values. The initial estimates of these hyperparameters should be refined during the repeated search procedure. In the rth repetition (r>1), assuming that is the optimal support chosen in the previous round and is the corresponding MAP estimate, we update the hyperparameters as follows:

4.3 Pilot Design

The performance of the proposed estimation scheme is mainly influenced by two issues: one is the approximation error which is equivalent to the data-induced ICI; the other is the reconstruction error. To eliminate the former, the pilots should be put together and far away from the data symbols. This consecutive pilot pattern enables to be extremely approximate to ՓP, thereby minimizing the approximation error. On the other hand, for the sparse signal reconstruction, the pilot symbols are suggested to be set at random places [4, 6]. However, in this case the pilots are very likely to be separated from each other and close to the data symbols.

For the trade-off between the accurate approximation of the system model and the stable reconstruction performance, we adopt the combination of the above two pilot patterns. Specifically, we divide the pilots into several groups and put these groups at random positions. Because each symbol mainly affects its adjacent ones, only the pilots at the borders of the groups are interfered by the unknown data symbols if the group size is sufficiently large. In this case, the data-induced ICI can be significantly reduced. Hence, compared with the random pilot pattern, the proposed combo pilot pattern has an improved robustness to the ICI. On the other hand, compared with the consecutive pilot pattern, the combo pilot pattern is expected to be more appropriate for the sparse signal reconstruction. The actual effects of the three pilot patterns are investigated in the next section.

 

5. Numerical Simulations

5.1 Simulation Settings

We assume that the OFDM system is divided into 1024 subcarriers, and the CP length ratio is 1/ 4 . We also simulate a doubly selective channel with ten paths, whose magnitudes α1,...,α10 are drawn from CN (0,1) . The delay-Doppler points (τs/T, fsNT) are uniformly chosen from the subset {0,1,...,255}×[-fmaxNT, fmaxNT], where the maximum normalized Doppler frequency fmaxNTvaries from 0.1 to 0.3. The additive channel noise is a zero-mean Gaussian process. The pilots are independently and randomly sampled from the QPSK alphabets, and the pilot ratio is 1 / 8.

Based on the approximate measurement model (34) with D = 2 , we reconstruct u using the modified BMP algorithm. The support search is repeated R = 3 times, and both the AMMSE and MAP estimates are calculated. The initial value of σ2 is roughly estimated using the LS estimate [1]. The initialization of the other hyperparameters is given by: Smax = 7 , P = Smax / 256 and = σ2/10. Moreover, we present the performances of the basis pursuit de-noising (BPDN) algorithm [19] and the orthogonal matching pursuit (OMP) algorithm [20] for comparison. Finally, the runtime of the estimation is also investigated.

5.2 Numerical Results

We first focus on the comparison between the proposed Bayesian approach and the conventional CS-based methods. The combo pilot pattern is adopted with the group size of 2D+1 . Fig. 2 depicts the mean square error (MSE) performance versus SNR when fmaxNT=0.2. Fig. 3 shows the MSE performance versus fmaxNT when SNR = 20 dB . On different conditions, the performances of the two modified-BMP estimators (AMMSE and MAP) are clearly superior to those of BPDN and OMP. The performance advantage demonstrates that the proposed Bayesian method effectively exploits the prior statistics and the group-sparsity of channel coefficients. Compared with the MAP estimate, the AMMSE estimate has an improved performance as expected. Moreover, the gap between the two becomes larger when low SNR or large Doppler spreading occurs, indicating that the AMMSE estimate is more robust to extreme channel conditions. This additional superiority of AMMSE is because that the repeated support search procedure enhances its performance in poor communication environments.

Fig. 2.Performance versus the SNR.

Fig. 3.Performance versus the Doppler frequency.

Next, Figs. 4 and 5 depict the performance of the modified BMP algorithm for different pilot patterns. Here, we compare the three previously mentioned pilot patterns: Pattern 1 is the consecutive pilot pattern; Pattern 2 is the random pilot pattern; Pattern 3 is the combo pilot pattern. The AMMSE estimate is adopted. It is observed that Pattern 3 outperforms the other two for different SNRs and Doppler frequencies. The good performance of this pattern demonstrates that it effectively reduces both the data-induced ICI and the reconstruction error. From Fig. 5, it is seen that the performance of Pattern 2 is close to that of Pattern 3 if the Doppler spreading is small. However, the former is more sensitive to fmaxNTs since it cannot address the severe ICI. Pattern 1 has a much worse performance than the other two pilot patterns do, indicating that it is not applicable for sparse recovery methods. In summary, the combo pilot pattern is very likely to be optimal for the proposed Bayesian estimation scheme over rapidly time-varying channels.

Fig. 4.Performances of different pilot patterns versus the SNR

Fig. 5.Performances of different pilot patterns versus the Doppler frequency

Finally, Figs. 6 and 7 show the mean runtime of these recovery methods. Here we still adopt the AMMSE estimate in the two BMP methods. As can be seen in the plots, the two BMP algorithms need shorter runtime than does the traditional CS-based method because of the efficient use of channel statistics. It is also observed that the modified BMP method is much faster than the original one, demonstrating that the proposed modifications effectively simplify the reconstruction steps.

Fig. 6.Mean runtime versus the SNR.

Fig. 7.Mean runtime versus the Doppler frequency.

 

6. Conclusion

We considered the application of the BMP algorithm to the channel estimation for rapidly time-varying OFDM systems. By exploiting the basis expansion function of the channel, we depicted the group-sparse structure and statistical characteristics of the doubly selective channel. Then, we proposed a BMP-based estimation scheme to reconstruct the CSI. In this method, we adopted both the MMSE and MAP estimates, providing a trade-off between the estimate performance and the computational complexity. Moreover, because the BMP algorithm cannot take advantage of the specific structure of the sparsity pattern, we made some changes to the support search procedure of BMP, which enable the algorithm to address the group-sparse channel coefficients in a more efficient manner. We showed that the modified BMP algorithm has a significantly reduced complexity compared with its original version. Furthermore, a pilot arrangement was designed to work in conjunction with the proposed estimation scheme. Experimental results illustrated the superiority of our method over conventional CS-based approaches.