Novel Adaptive Distributed Compressed Sensing Algorithm for Estimating Channels in Doubly-Selective Fading OFDM Systems

Abstract

Doubly-selective (DS) fading channel is often occurred in many orthogonal frequency division multiplexing (OFDM) communication systems, such as high-speed rail communication systems and underwater acoustic (UWA) wireless networks. It is challenging to provide an accurate and fast estimation over the doubly-selective channel, due to the strong Doppler shift. This paper addresses the doubly selective channel estimation problem based on complex exponential basis expansion model (CE-BEM) in OFDM systems from the perspective of distributed compressive sensing (DCS). We propose a novel DCS-based improved sparsity adaptive matching pursuit (DCS-IMSAMP) algorithm. The advantage of the proposed algorithm is that it can exploit the joint channel sparsity information using dynamic threshold, variable step size and tailoring mechanism. Simulation results show that the proposed algorithm achieves 5dB performance gain with faster operation speed, in comparison with traditional DCS-based sparsity adaptive matching pursuit (DCS-SAMP) algorithm.

1. Introduction

Nowadays, wireless communications are required to provide high-speed and reliable data transmission in high-speed mobile environments, which unfortunately leads to doubly-selective (DS) fading in orthogonal frequency division multiplexing (OFDM) systems, such as high-speed rail communication systems and underwater acoustic (UWA) wireless networks. Due to the coherent receive under the condition of fast fading, it is urgent to develop fast channel estimate techniques over DS fading OFDM channels [1].

There are many existing works [2], [3] verifying the sparsity of DS channel. And the authors in [3] also illustrate that DS channel can be approximated by basis expansion model (BEM). In [4], a novel channel estimator based on BEM is employed in IEEE 802.16e systems under high-speed mobile environment, which significantly reduces the estimated channel coefficients. However, the sparsity of channel is not exploited to improve the estimation performance. And also combined with BEM, a compressive sensing (CS) algorithm is proposed in [5] to estimate DS channel which performs better but requires much less pilots than traditional estimations, such as least squares (LS) method. However, the CS-based method recovers each sparse signal individually, which is inefficient in the reconstruction of multiple correlated signals.

Distributed compressive sensing (DCS) [6] theory has been developed to jointly reconstruct the set of sparse signals by exploiting their joint sparsity. A set of sparse signals exhibit joint sparsity, which is characterized by that the non-zero elements in each signal vector have the same location indices but different values. In [7], the DCS simultaneous orthogonal matching pursuit (DCS-SOMP) algorithm is proposed and employed to estimate the DS channel based on complex exponential-BEM (CE-BEM), which is superior to the CS method. Furthermore, the work [8] improves the SOMP estimation by adopting exponential smoothing, which gets better performance under the UWA communication. However, these estimation methods need the channel sparseness as a prior information, which is often not available in the practical communication system. In order to overcome this drawback, the DCS-sparsity adaptive matching pursuit (DCS-SAMP) algorithm [9] has been proposed, which can reconstruct signals without any prior information about the sparsity of signals. But there are still deficiencies in its accuracy of reconstruction and speed of operation.

In this paper, we propose an improved DCS algorithm, which is called DCS-improved sparsity adaptive matching pursuit (DCS-IMSAMP) for estimating DS channel without the knowledge of channel sparsity. In the proposed algorithm, the dynamic threshold which is motived by stagewise weak orthogonal matching pursuit (SWOMP) [10] is employed in prelim test to enhance the accuracy of element selections per iteration. Furthermore, the variable step size and tailoring mechanism are adopted to improve operation speed.

Simulations are conducted to show that the proposed channel estimation method outperforms the scheme of DCS-SAMP about 5 dB and costs less computation time under the condition of unknown sparsity.

The reminder of this paper is organized as follows. In Section 2, the DS channel estimation model in OFDM systems and the BEM theory are introduced. The DCS-SAMP and DCS-IMSAMP algorithm are depicted in Section 3. In Section 4, simulation results are displayed to illustrate the superior performance of our proposed DCS-IMSAMP algorithm, as compared with existing DCS-SAMP algorithm. Finally, we summarize our work in Section 5.

Notations: For a matrix \(\mathbf{\Phi}, \mathbf{\Phi}^{-1}, \mathbf{\Phi}^{H} \text { and } \mathbf{\Phi}^{+}\) denote inverse, conjugate transpose and pseudo inverse, respectively. \(\|\mathbf{\Phi}\|_{r o w-\infty}\) denotes taking ∞ -norm by row of matrix Φ . Matrix \(\mathbf{\Phi}_{C}\) is the submatrix of Φ with column indices in the set C and all rows. The function \(\operatorname{Max}\{\varphi, U\}\) acts to return U indices corresponding to the largest values of the vector φ .

2. OFDM System Over a Doubly Selective Channel Model Based on BEM

Considering an OFDM system with N subcarriers and L channel taps. Vector \(\mathbf{h}_{l}=(h[0, l], \ldots, h[N-1, l])^{T} \in C^{N \times 1}\) denotes the time variation of the l-th \((l \in[0, L-1])\) channel tap in a sampling period \(T_s\) , where \(h[n, l]\) denotes the l-th discrete channel impulse response at time \(n(n \in[0, N-1])\) as

\(h[n, l]=\sum_{i=0}^{I-1} \eta_{i} e^{j 2 \pi v_{i} n T_{\text {sinc }}}\left(l-\frac{\tau_{i}}{T_{s}}\right),\)       (1)

where \(\eta_{i}, v_{i} \text { and } \tau_{i}\) respectively denote the complex attenuation factor, Doppler frequency, and time delay of i-th path. The BEM model describes the time variation of DS channel within a certain time by linear combination of finite basis functions [11], which can efficiently reduce the channel coefficients to be estimated. Using the BEM, the l-th channel tap \(h_l\) can be written as

\(\mathbf{h}_{l}=\left(\mathbf{b}_{0}, \mathbf{b}_{1}, \ldots, \mathbf{b}_{Q-1}\right)\left(\begin{array}{c} z[0, l] \\ z[1, l] \\ \vdots \\ z[Q-1, l] \end{array}\right)+\xi_{l},\)       (2)

where \(\mathbf{b}_{q}=(b[0, q], \ldots, b[N-1, q])^{T}\) is the q-th (q∈[0,Q-1]) BEM function, z[q,l] is the BEM coefficients of corresponding basis function, \(\xi_{l}=(\xi[0, l], \ldots, \xi[N-1, l])^{T}\) represents model error which will be ignored in the following discussion, and \(Q(Q<<N)\) denotes the BEM order which is generally defined as an odd number. According to (2), it is clearly that BEM method can reduce the total number of undetermined coefficients from NL to QL. We utilize the CE-BEM to approximate the DS channel, where basis function is written as \(b_{q}(n)=e^{j \omega_{q} n}, \quad \omega_{q}=\frac{2 \pi}{N}\left(q-\frac{Q-1}{2}\right)\).

The DS channel is sparse in the delay domain, that is, only K tap coefficients have larger values and other tap coefficients can be approximately zero, due to the fast decay of sinc(·) function [12]. Assuming \(h(n, l)=0, l \notin \mathbf{\Gamma}\), where \(\Gamma=\left\{f_{1}, f_{2}, \ldots f_{K}\right\}\) denotes the set of positions for the K larger channel taps. Basing on this, \(z[q, l]=0, l \notin \mathbf{\Gamma}\) can be derived by equation (2). Then define \(\mathbf{z}_{q}=(z[0,0], \ldots, z[0, L-1])^{T} \in C^{L x 1}\) as the coefficient vector of the q-th CE-BEM basis function. Obviously, the column vectors \(\left\{\begin{array}{l} z \\ q \end{array}\right\} \begin{array}{l} Q-1 \\ q=0 \end{array}\) have the same non-zero element locations but different non-zero values, which are jointly sparse.

Hence, the received signal \(\mathbf{Y} \in C^{N \times 1}\) can be expressed as

\(\begin{aligned} \mathbf{Y}=& \mathbf{H}_{F} \mathbf{X}+\mathbf{W} \\ &=\left(\sum_{q=0}^{Q-1} \mathbf{F}_{N} \operatorname{diag}\left(\mathbf{b}_{q}\right) \mathbf{F}_{N}^{H} \mathbf{Z}_{q}\right) \mathbf{X}+\mathbf{W} \end{aligned}\)       (3)

where \(\mathbf{H}_{F}\) is the N×N channel matrix in frequency domain, X is the N×1 transmit signal vector, \(\mathbf{F}_{N}\) is the L×N normalized IFFT matrix with (l,m)-th entry expressed as \(\left[\mathbf{F}_{N}\right]_{l, m}=\exp (-j 2 \pi l m / N)\), \(\mathbf{Z}_{q}=\operatorname{diag}\left(\sqrt{N} \mathbf{F}_{N}\left(\mathbf{z}_{q}^{T}, 0_{1 \times(N-L)}\right)^{T}\right)\) is the N×N diagonal matrix.

The sparse pilot pattern, which consists of non-zero pilots and guard pilots, is employed in the OFDM system. The total number of pilot subcarriers is \((2 Q-1) T(T<<L)\) and let P denotes the set of indices for pilots. Pilot subcarriers are divided into T clusters, including one non-zero pilot and (2Q-2) guard pilots which can effectively reduce inter-subcarrier interference in high-speed mobile scenarios. The constant amplitude zero auto correlation sequences are used as the transmitted data on the non-zero pilot subcarriers in this paper, which is denoted as

\([\mathbf{X}]_{r_{r}}=\mathbf{P}_{r a l}=\left([\mathbf{X}]_{p_{0}}, \dots,[\mathbf{X}]_{p_{r-1}}\right)^{T},\)       (4)

where\(P_{v}=\left\{p_{0}, \ldots, p_{T-1}\right\} \quad\left(0 \leq p_{0}<\ldots p_{T-1} \leq N-1\right)\) is the set of non-zero pilot indices. And the data of guard pilot and corresponding index set are respectively expressed as \([\mathbf{X}]_{P_{g}}=\mathbf{0}\) and

\(P_{g}=\left\{p_{r}-Q+1\right\} \cup \ldots \cup\left\{p_{r}-1\right\} \cup \ldots \cup\left\{p_{r}+Q-1\right\},\)       (5)

In order to conform to DCS model, all pilot indices in one OFDM symbol are re-divided into Q groups, which is written as        

\(\left\{P_{0}, \ldots, P_{Q-1}\right\}=\left\{P_{v}-\frac{Q-1}{2}, \ldots, P_{v}, \ldots, P_{v}+\frac{Q-1}{2}\right\}\)       (6)

And the arrangement of pilot subcarriers with Q=3 is presented in Fig. 1. There is a non-zero pilot marked as a solid circle in the middle of each cluster and there are two zero-valued guard pilots marked as squares on each side of this non-zero pilot. While, other circles represent as data subcarriers which carry the transmit data.

 

Fig. 1. Arrangement of pilot subcarriers (Q=3)

Consequently, the received pilot signals are expressed as

,       (7)

where \(\begin{equation}\left[\mathbf{F}_{N}^{\prime}\right]_{P_{\frac{Q-1}{2}}} \in C^{T \times L}\end{equation}\) denotes the submatrix of \(\begin{equation}\mathbf{F}_{N}^{\prime}\end{equation}\) with row indices \(\begin{equation}P_{\frac{Q-1}{2}}\end{equation}\) and columns from 1 to L, \(\begin{equation}\left\{\mathbf{W}_{q}\right\}_{q=0}^{Q-1} \in C^{T \times I}\end{equation}\) includes model errors and noises. We further define \(\begin{equation}\mathbf{\Phi}=\operatorname{diag}\left(\mathbf{P}_{\text {ral}}\right)\left[\mathbf{F}_{N}\right]_{\frac{p-1}{2}} \in C^{T \times L}\end{equation}\), \(\begin{equation}\mathbf{Z}=\left[\begin{array}{lll} \mathbf{z}_{0} & \ldots \mathbf{z}_{\frac{Q-1}{2}} & \ldots \mathbf{z}_{Q-1} \end{array}\right] \in C^{L \times Q}\end{equation}\), received pilot signals \(\begin{equation}\mathbf{Y}_{P}=\left[[\mathbf{Y}]_{P_{0}} \ldots[\mathbf{Y}]_{P_{\frac{Q-1}{2}}} \ldots[\mathbf{Y}]_{P_{Q-1}}\right] \in C^{T \times Q}\end{equation}\)

, and \(\begin{equation}\mathbf{W}=[\mathbf{W}_{0} \ldots \mathbf{W}_{\frac{Q-1}{2}} \ldots \mathbf{W}_{Q-1}] \in C^{T \times Q} \end{equation}\) .

Then, according to (7), the received pilot matrix can be written as      

\(\begin{equation}\mathbf{Y}_{P}=\mathbf{\Phi Z}+\mathbf{W}\end{equation}\)       (8)

3. Channel Estimation Algorithm

DCS theory provides different joint sparsity models according to different actual scenarios, including the common spare supports model (JSM-2) [13]. In this model, all sparse signals \(\begin{equation}\boldsymbol{\theta}_{q} \in C^{L \times 1}, \forall q \in\{0,1, \ldots, Q-1\}\end{equation}\) have K nonzero entries, which own different values but are in the same positions in each θ_q . The objective of DCS reconstruction algorithm is to jointly reconstruct the \(\begin{equation}\left\{\boldsymbol{\theta}_{q}\right\}_{q=0}^{Q-1}\end{equation}\) based on the obtained observation matrix \(\begin{equation}\left\{\mathbf{s}_{q}\right\}_{q=0}^{Q-1}\end{equation}\) and a known measurement matrix \(\begin{equation}\mathbf{\Psi} \in C^{T \times L}\end{equation}\) , which is defined as      

\(\begin{equation}\mathbf{s}_{q}=\mathbf{\Psi} \boldsymbol{\theta}_{q}, \forall q \in\{0,1, \ldots, Q-1\}\end{equation},\)       (9)

Accordingly, doubly selective OFDM channel model based on CE-BEM can be regarded as a JSM-2 with additive noise, which can be written as        

\(\begin{equation}\left[\mathbf{Y}_{P}\right]_{q}=\mathbf{\Phi} \mathbf{z}_{q}+\mathbf{W}_{q}, \forall q \in\{0,1, \ldots, Q-1\}\end{equation},\)       (10)

So, we can apply DCS reconstruction algorithm to complete the doubly selective OFDM channel estimation, which clearly outperforms the scheme based on CS theory according to [7], [14]. The DCS-SOMP algorithm has been used to estimate the jointly sparse channel in [7], which achieves much better reconstruction performance than CS-based method. However, the DCS-SOMP-based estimation scheme is hardly employed in practical cases, since it is very difficult to obtain the accurate sparsity of channel as a prior condition. Thus, it is necessary to study the channel estimation method when the sparsity of channel is unknown. The DCS-SAMP has been proposed to solve this problem, while it still has flaws in accuracy and operation speed.

3.1 DCS-SAMP Algorithm

DCS-SAMP can reconstruct jointly sparse signals without the prior knowledge of the signal’s sparsity. It progressively estimates the true sparsity and the support set of the recovery signals. The DCS-SAMP-based channel estimation is as follows:

Input: measurement matrix \(\begin{equation}\mathbf{\Phi} \in C^{T \times L}\end{equation}\) , received pilot matrix \(\begin{equation}\mathbf{Y}_{P} \in C^{T \times Q}\end{equation}\) , initial step S=1;

Output: the estimated CE-BEM coefficient matrix \(\begin{equation}\hat{\mathbf{Z}} \in C^{L \times Q}\end{equation}\) ;

 1) Initialization:

 \(\left.\mathbf{r}_{0}=\mathbf{Y}_{p} \quad \text { \{initial residue }\right\}\\ \Lambda_{0}=\varnothing \quad\{\text { the estimated final list }\} \\ U=S \quad\{\text { size of the support set in the first stage }\}\\ k=1 \quad \text { \{ iteration index }\}\\ j=1 \quad \{\text { stage index }\}\)

 2) Iterate:

 a) \(\left.J_{k}=\operatorname{Max}\left\{\left\|\Phi^{H} \mathbf{r}_{k-1}\right\|_{\text {Irow-oo }}, U\right\} \quad \text { \{Prelim Test }\right\}\)

 b) \(\left.C_{k}=\Lambda_{k-1} \cup J_{k} \quad \text { \{Update Candidate List }\right\}\)

 c) \(\left.F=\operatorname{Max}\left\{\left\|\mathbf{\Phi}_{C}^{+} Y_{P}\right\|_{\text {row-oo }}, U\right\} \quad \text { \{Final Test }\right\}\)

 d) \( \mathbf{r}_{n e w}=\mathbf{Y}-\mathbf{\Phi}_{F}\left(\mathbf{\Phi}_{F}^{H} \mathbf{\Phi}_{F}\right)^{-1} \mathbf{\Phi}_{F}^{H} \mathbf{Y}_{P} \quad\{\text { Calculate Residual }\}\)

 e) If halting condition is satisfied, then

 stop the iteration and go to step 3);

 Else if \(\left\|\mathbf{r}_{\text {new}}\right\|_{2} \geq\left\|\mathbf{r}_{k-1}\right\|_{2}\) then

 \( j=j+1, \quad U=j \times S;\)

 Else

 \( \mathbf{r}_{k}=\mathbf{r}_{n e w}, \quad \Lambda_{k}=F, \quad k=k+1;\)

 End if

 Until the stopping condition holds

 3) Calculate:  \(\overline{\mathbf{Z}}_{\Lambda_{k}}=\left(\mathbf{\Phi}_{\Lambda_{k}}^{H} \mathbf{\Phi}_{\Lambda_{k}}\right)^{-1} \mathbf{\Phi}_{\Lambda_{k}}^{H} \mathbf{Y}_{P}\)

 4) Output: The estimated coefficient matrix \(\hat{Z}\) which has non-zero values at the position of set Λ_k , the corresponding value is \(\overline{\mathbf{Z}}_{\Lambda_{k}}\) and the rest are zero. Here, set C and Λ_k are defined respectively as candidate list and final list. After recovering coefficient matrix \(\hat{\mathbf{Z}}=\left[\begin{array}{lll} \hat{\mathbf{z}}_{0} & \ldots & \hat{\mathbf{z}}_{\frac{Q-1}{2}} & \ldots \hat{\mathbf{z}}_{Q-1} \end{array}\right]\)

\(\{\hat{\mathbf{h}}\}_{l=0}^{L}\) can be obtained according to equation (2).

However, the algorithm has deficiencies in accuracy of the prelim test and choice of its initial step. In the prelim test, U elements are selected in each iteration. But, with the increment of the stage index j, U will also increase, which results in the selection of too many elements. Since candidate list contains extra elements, the accuracy of the selection is reduced significantly, and more time is wasted on the calculation of the misjudgment elements. Furthermore, the value of the initial step is very important for the algorithm's estimation precision and running time. In order to ensure the accuracy of the algorithm, generally, initial step size S is 1, which results in an increase of the running time. Thus, the DCS-SAMP obviously needs to be improved.

3.2 Our Proposed DCS-IMSAMP Algorithm

In this paper, a new algorithm named DCS-improved sparsity adaptive matching pursuit (DCS-IMSAMP) is proposed which improves the operation speed and accuracy without the prior information about the sparsity of channel, as compared with the DCS-SAMP algorithm. The proposed algorithm utilizes the joint sparsity of signals and adopts the dynamic threshold idea in the SWOMP algorithm as preliminary rule to improve the accuracy of the algorithm, and improves the speed by choosing variable step size and tailoring mechanism. The main steps of the algorithm are as follows:

Input: measurement matrix \(\mathbf{\Phi} \in C^{T \times L}\) , received pilot matrix \(\mathbf{Y}_{P} \in C^{T \times Q}\) , initial step S;

Output: the estimated CE-BEM coefficient matrix \(\hat{\mathbf{Z}} \in C^{L \times Q}\) ;

 1) Initialization:

 \(\left.\mathbf{r}_{0}=\mathbf{Y}_{P} \quad \text { \{initial residue }\right\}\\ \Lambda_{0}=\varnothing \quad\{\text { the estimated final list }\}\\U=S \quad\{\text { size of the support set in the first stage }\}\\ k=1 \quad \text { \{iteration index }\}\\ I=0 \quad\{\text { stage index }\}\)

 2) Iterate:

 a) \(J_{k}=\left\{j:\left|\mathbf{u}_{k}(j)\right|>\alpha \cdot u_{\max }\right\}, where \mathbf{u}_{k}=a b s\left[\left|\Phi^{H} r_{k-1}\right|_{r o u s-\infty}\right] \quad\{\text { Prelim Test }\}\)

 b) \(C_{k}=\Lambda_{k-1} \cup J_{k} \quad\{\text { Update Candidate List }\}\)

 c) \( If \left.\operatorname{size}\left(C_{k}\right)>\mu \cdot T \text { then } I=1 \text { End if } \quad \text { (Determine Stage Index }\right\}\)

 d) \(If \operatorname{size}\left(C_{k}\right) \geq U then F=\operatorname{Max}\left(\left|\Phi_{C_{k}}^{+} \mathbf{Y}_{P_{P o u s-\infty}}\right|_{\mathbf{r}}, U\right) \quad\{\text { Final Test }\}\)

 Else \(F=C_{k};\)

 End if

 e) \( \mathbf{r}_{n e w}=\mathbf{Y}_{p}-\mathbf{\Phi}_{F}\left(\mathbf{\Phi}_{F}^{H} \mathbf{\Phi}_{F}\right)^{-1} \mathbf{\Phi}_{F}^{H} \mathbf{Y}_{P} \quad\{\text { Calculate Residual }\}\)

 f) if \(\left\|\mathbf{r}_{\text {new}}\right\|_{2}<\varepsilon\) then

 stop the iteration and go to step 3);

 Else if \(\left\|\mathbf{r}_{\text {new}}\right\|_{2} \geq\left\|\mathbf{r}_{k-1}\right\|_{2}\) and \(I=0\) then

 \(U=U+2 \times S, \Lambda_{k}=\Lambda_{k-1}, \quad \mathbf{r}_{k}=\mathbf{r}_{k-1};\)

 Else if \(\left\|\mathbf{r}_{\text {rew }}\right\|_{2} \geq\left\|\mathbf{r}_{k-1}\right\|_{2}\) and \(I=1\) then

 \(U=U+S, \Lambda_{k}=\Lambda_{k-1}, \mathbf{r}_{k}=\mathbf{r}_{k-1};\)

 Else

 \(\mathbf{r}_{k}=\mathbf{r}_{n e w}, \quad \Lambda_{k}=F;\)

 End if

 g) k=k+1 {Update Iteration Index}

 Until the stopping condition holds

 3) Calculate: \(\overline{\mathbf{Z}}_{\Lambda_{k}}=\left(\mathbf{\Phi}_{\Lambda_{k}}^{H} \mathbf{\Phi}_{\Lambda_{k}}\right)^{-1} \mathbf{\Phi}_{\Lambda_{k}}^{H} \mathbf{Y}_{P}\)

 4) Output: The estimated coefficient matrix \(\hat{\mathbf{Z}}\) which has non-zero values at the position of set \(\Lambda_{k}\) , the corresponding value is \(\overline{\mathbf{Z}}_{\Lambda_{k}}\) and the rest are zero.

Consequently, we can obtain \(\{\hat{\mathbf{h}}\}_{l=0}^{L}\) based on \(\hat{\mathbf{Z}}=\left[\begin{array}{lll} \hat{\mathbf{Z}}_{0} & \ldots \hat{\mathbf{z}}_{\frac{Q-1}{2}} & \cdots \hat{\mathbf{z}}_{Q-1} \end{array}\right]\) . In the proposed algorithm, there are three mainly differences from DCS-SAMP algorithm. Firstly, the preliminary threshold is changed, which is based on the maximum of current \(u_k\) , where a∈[0.5,1] is a preset threshold parameter and \(u_{max}\) is the maximum element in the vector \(u_k\) . Instead of selecting U elements in preliminary test, we propose to merely select elements satisfying new condition in each iteration, which is inspired by the SWOMP. Using this idea, misjudgment elements are reduced.

Secondly, unlike the DCS-SAMP algorithm based merely on backtracking, our proposed DCS-IMSAMP algorithm combines the tailoring mechanism in the final test to remove misjudgment coordinates. After updating the candidate list C_k , in step 2(d), the columns of the measurement matrix corresponding to the element of the updated candidate set are taken out to form a sub-matrix \(\Phi_{C_{k}}\) , which then forms a least square solution with the matrix \(\mathbf{Y}_{P}\) . And using function \(\operatorname{Max}\left(\left|\mathbf{\Phi}_{C_{k}}^{+} \mathbf{Y}_{P}\right|_{\text {rows}-\infty}, U\right)\) to find the U coordinates of the largest values to remove other incorrect coordinates, if the number of indices in C_k is larger than U.

Lastly, the variable step size is adopted to approach the true sparsity in the proposed algorithm, which obtains the better adaptability than the fix step size used in DCS-SAMP. In step 2(f), the algorithm updates the step size basing on the current residual and stage index I which is determined in step 2(c). And µ<1/4, according to the 1/4 practical rule [15].

When I = 0, it indicates that the size of candidate list is small, so the next iteration takes a big step to effectively expand the support set; I = 1, which indicates that the size of candidate list is large. And in order to prevent over-estimation of sparsity, a small step is used to filter the atoms and approximate the sparsity. Therefore, the initial step size S will be more flexible.

4. Simulation Results

In this section, three different algorithms, such as our proposed DCS-IMSAMP algorithm, DCS-SOMP algorithm and DCS-SAMP algorithm, are employed to complete the DS OFDM channel estimation under the condition that the sparsity of channel is unknown. Monte Carlo simulations are conducted to compare the performances of the DS channel estimation employing above three algorithms. To evaluate the estimation performance, the normalize mean square error (NMSE) is used, which is defined as

\(\operatorname{NMSE}(d B)=10 \log _{10}\left(\frac{\mathrm{E}\left(\|\mathbf{h}-\hat{\mathbf{h}}\|_{2}^{2}\right)}{\mathrm{E}\left(\|\mathbf{h}\|_{2}^{2}\right)}\right),\)       (11)

where h denotes the true CIR and \(\hat{h}\) is the estimated CIR. The parameters of CE-BEM model and OFDM system are shown in Table 1.

Table 1. Parameters of simulation 

 

4.1 Research on the Threshold Parameter a for DCS-IMSAMP

In this experiment, we explore the optimal selection of a in the channel estimation by simulations, since the parameter a takes important role in the preliminary test. We examine different values of a, such as 0.5, 0.55, 0.65, 0.77, 0.85. Fig. 2 displays the NMSE performances with different a, as can be showed, when  a ∈ [0.50,0.65], the performances of channel estimation are superior to those with other values of a.

To further compare the operation speed of the schemes for a between 0.5 and 0.85, Table 2 presents the average running time for a single channel estimation and their NMSE at SNR=18dB. It can be found that the running time increases, as the value of a increases.

However, there is only a tiny difference in running time, when a ∈ [0.50,0.65]. Therefore, combined with results in Fig. 2, the parameter a is fixed at 0.55 in the following stimulations.

 

Fig. 2. NMSE performances with different a under unknown sparsity.

Table 2. Average running times for different threshold parameters a

 

4.2 Research on the Initial Step Size S for DCS-IMSAMP

As above discussion in Section 3.2, initial step size S can be chosen more flexibly in the DCS-IMSAMP algorithm. However, in the DCS-SAMP algorithm, the S should not be too large, generally taken as 1, to ensure good channel estimation performance. Fig. 3 shows the comparison of the NMSE performance between DCS-IMSAMP with different initial steps and DCS-SAMP with S = 1. Obviously, for DCS-IMSAMP, the smaller initial step size is, the lower NMSE gets. Particularly, the performances of channel estimation are similar when initial steps are taken from 1 to 4, especially with SNR ∈ [25,35]. Meanwhile, the performance of DCS-IMSAMP outperforms the DCS-SAMP with S = 1, even with large initial step. In the Section 4.4, it will be verified that the algorithm takes less running time when initial step is chosen larger value.

 

Fig. 3. NMSE performance between DCS-IMSAMP with S = 1,4,8,12 and DCS-SAMP with S = 1 under unknown sparsity

4.3 NMSE Comparison Between DCS-SOMP, DCS-SAMP and DCS-IMSAMP

To illustrate the better performance of our proposed scheme, the DCS-SOMP-based, DCS-SAMP-based and DCS-IMSAMP-based channel estimation methods are conducted respectively under the same channel states and unknown sparsity of channel. For the sake of the fairness of comparison, the stopping iteration condition of the DCS-SOMP algorithm is changed from the maximum number of iterations to a threshold parameter \(\varepsilon=10^{-6}\). And the initial steps of DCS-SAMP and DCS-IMSAMP are both taken as 1.

In Fig. 4(a), we make the NMSE comparison between three schemes in the condition of normalized Doppler frequencies  \(v_{\max }=0.065\). It is clearly shown that the proposed algorithm achieves the higher accuracy of channel estimation than other algorithms. For example, when the SNR is in the range of 5 to 20, the NMSE of DCS-IMSAMP algorithm is at least 5 dB lower than that of the DCS-SAMP scheme. On the other hand, the DCS-SOMP-based scheme has the worst performance because it approximates the actual sparsity of channel only by judging the size of the residual.

 

Fig. 4. NMSE performance between DCS-SOMP, DCS-SAMP, DCS-IMSAMP with \(v_{\max }=0.065\), 0.093 under unknown sparsity

In order to further demonstrate the superiority of the proposed algorithm, Fig. 4(b) plots the performance of these schemes under \(v_{\max }=0.093\). Again, the simulation results show that the performance of DCS-IMSAP far exceeds that of others. The performance of channel estimation in all algorithms decline, mainly because the CE-BEM modeling error increases with the increasement of Doppler frequency.

4.4 Computational Complexity Comparison of Different Channel Estimation Methods

In this subsection, for each scheme, we perform 800 simulations and calculate the average time taken for a single channel estimate at different SNRs to compare their computational complexity. The results of different algorithms are shown in Fig. 5.

 

Fig. 5. Comparison of running time between DCS-SOMP, DCS-SAMP and DCS-IMSAMP (S=1,4,10)

It can be observed that channel estimation under unknown sparsity by DCS-SOMP algorithm spends the most time, due to the change of its stopping condition. Clearly, the running time by DCS-IMSAMP is less than that of DCS-SAMP, regardless of whether its initial step takes a larger value or a smaller value. This is because that the adoption of variable step size and tailoring mechanism can quickly approximate the true sparsity and save running time effectively. And for DCS-IMSAMP, the running time gets shorter as the step size increases. Thus, combined with Fig. 4(a), the initial step in DCS-IMSAMP can be appropriately chosen to be larger, such as S=4, so as to reduce the running time and have less influence on the accuracy of channel estimation.

5. Conclusion

In this paper, the DCS-IMSAMP algorithm was proposed to complete the channel estimation of DS channel in OFDM systems. Compared with other DCS algorithms, such as DCS-SOMP, DCS-SAMP, the advantages of DCS-IMSAP are its higher accuracy of channel estimation and faster operation speed in the condition of unknown sparsity. Meanwhile, the proposed algorithm can also be applied in other scenarios, like detection of massive MIMO, image recovery and so on.