Spatially Correlated Sparse MIMO Channel Path Delay Estimation in Scattering Environments Based on Signal Subspace Tracking

A parametric scheme for spatially correlated sparse multiple-input multiple-output (MIMO) channel path delay estimation in scattering environments is presented in this paper. In MIMO outdoor communication scenarios, channel impulse responses (CIRs) of different transmit–receive antenna pairs are often supposed to be sparse due to a few significant scatterers, and share a common sparse pattern, such that path delays are assumed to be equal for every transmit–receive antenna pair. In some existing works, an exact common support condition is exploited, where the path delays are considered equal for every transmit–receive antenna pair, meanwhile ignoring the influence of scattering. A more realistic channel model is proposed in this paper, where due to scatterers in the environment, the received signals are modeled as clusters of multi-rays around a nominal or mean time delay at different antenna elements, resulting in a non-strictly exact common support phenomenon. A method for estimating the channel mean path delays is then derived based on the subspace approach, and the tracking of the effective dimension of the signal subspace that changes due to the wireless environment. The proposed method shows an improved channel mean path delays estimation performance in comparison with the conventional estimation methods.


Introduction
Multiple-input multiple-output (MIMO) technology has become an active research topic during the last decade due to its capability for achieving the high transmission rates required by an increasing number of data-demanding applications. MIMO technology provides a plenty of benefits that allow dealing with the challenges posed by the impairments in the wireless channel, especially multipath fading [1]. It provides several important performance gains such as antenna gain, diversity gain, and multiplexing gain. The benefits of MIMO are achieved through the exploitation of the spatial dimension across multiple antennas at the transmitter and receiver, in addition to the time and frequency dimensions already exploited in the conventional single-input single-output (SISO) systems [2]. However, a performant "ideal" MIMO communication system would require an exact knowledge of the MIMO channel or channel state information (CSI). MIMO channel parameters estimation is required for the equalization at the receiver side and precoding at the transmitter side. Hence, channel parameter estimation is a crucial part in MIMO communication. Channel estimation methods are divided into nonparametric and parametric categories. Nonparametric methods are unconstrained methods that exploit no or few assumptions about the channel model, such that they have no associated parameters to rely on, resulting in an increase in the dimensions of the estimation problem. On the other hand, parametric methods are associated with parametric models that are defined by a finite number of parameters with clear physical meanings related to the signal propagation through the channel. Among many other characteristics, some parametric methods exploit the sparsity of wireless channels [3][4][5]. Parametric methods, if wisely used, can achieve robust estimation and decrease the dimensions of the estimation problem.
An accurate MIMO channel model is important to exploit the advantages provided by MIMO systems. This should be done in a realistic way that accurately describes the propagation channel. Recently, the issue of sparsity of the wireless channel has received a lot of attention. It is a fact that, in many scenarios, wireless channels can be characterized by few significant paths due to a few dominant scatterers distributed in the environment, causing a sparse channel impulse response (CIR) [6][7][8][9][10]. This sparsity property has been shown to be reasonable through physical investigations that have been reported for some physical channels [11]. In MIMO outdoor communication scenarios, the CIRs of different transmit-receive antenna pairs share a common support. As the transmit and receive array dimensions are relatively small compared to the long transmission distance, the signals received by the different closely located antennas for each path are spatially correlated [9,12]. In some existing works [9,13], an exact common support condition is considered such that different transmit-receive antenna pairs share exactly the same path delays, meanwhile ignoring the effect of scattering due to the environment. Clustering due to scatterers is an important property that characterizes a wireless channel, such that the channel is modeled as clusters of multi-rays [14][15][16][17]. In this work, in order to match a more realistic channel model, a sparse clustered channel model is considered. The channel is assumed to be sparse in the sense that it contains few significant paths where, due to the physical properties of scatterers, each significant path is modeled as a cluster of multi-rays around a nominal or mean time delay [18][19][20]. This is what we refer to as "delay spreading". Considering the long distance between the transmit and receive antenna arrays in outdoor scenarios, the above mentioned delay spreading can be assumed to be relatively small. It follows that the multipath clusters can be basically described by their mean delays where it seems neither meaningful nor possible to estimate the time delay of each ray constituting the clusters.
In this paper, a new derivation of the channel model based on higher-order Taylor expansion is proposed. It is worth noting that some previous derivations of the channel model with scattering have been proposed in other contexts (direction of arrival estimation/localization, joint angle and delay estimation) but were limited to the first-order Taylor expansion [18][19][20]. A subspace based method that exploits the derived parametric channel model is then proposed to estimate the mean path delays. The minimum description length (MDL) criterion [21,22] is used to track the signal subspace in order to estimate its effective dimension that is changing due to the wireless environment. Taking into consideration this phenomenon, the proposed method provides an improved channel mean path delays estimation performance in comparison to the conventional subspace based methods which do not take the scattering in the environment into account.

Channel Model
For an N × M MIMO system, considering the models presented in [18][19][20] and focusing on the channel path delays, the CIR h (n,m) (t) between the nth transmit antenna and mth receive antenna can be modeled as: where δ(.) is the Dirac function; L is the number of propagation paths (clusters), with P l multi-rays for each path (cluster); (t l + τ (n,m) l p ) and α (n,m) l p denote the path delays and path gains of path l, respectively, between the nth transmit antenna and mth receive antenna, where t l is the mean delay associated to path l shared by all the transmit-receive antenna pairs, and τ (n,m) l p is the delay deviation associated to the pth contributing ray of path l between the nth transmit antenna and mth receive antenna.

System Model
In the considered N × M MIMO system, each antenna at the emitter side is transmitting pilot symbols on a different carrier in order to identify the channels associated with different transmit-receive antenna pairs. Consider a known pulse shape g(t) transmitted at the nth transmit antenna with a constant rate 1/T through a medium consisting of L paths. For the mth receive antenna, the received baseband signal is: For notational convenience, the additive noise is not represented in the above equation.
Applying the discrete Fourier transform (DFT), the Fourier coefficients of the received signal are given by: for k = −K/2 + 1, . . . , K/2, where K is the considered number of Fourier coefficients and G[k] is the DFT of pulse g(t). Note that K is an even integer.

Scattering Model with Non-Strictly Exact Common Support
Let v k (t) = e −j 2π T kt ; Equation (3) can be rewritten as: As the deviations τ (n,m) l p are considered small, the U th -order Taylor expansion of Equation (4) gives: where v (u) ) is the remaining term in the Taylor approximation, which is considered small, such that its impact in the approximation can be neglected.
Equation (5) can be written in the following form: where a (n,m) Equation (6) is written in the following matrix form: where It is worth noting that matrix V is the same for all the transmit-receive antenna pairs, as the mean delays of each path are assumed to be the same for all the transmit-receive antenna pairs.
The matrix of Fourier coefficient vectors at antenna m due to the N transmit antennas is given by: which can be written as: where A (m) = [a (1,m) . . . a (N,m) ].
Arranging the Fourier coefficient matrices of different receive antennas as follows: X can be written as: where Multiplying Equation (13) by the inverse of G, we have:

Second Order Statistical Analysis
The covariance matrix R Y of the observation Fourier coefficient vectors Y could be estimated as: For the noise-free case: where the covariance matrix R a is calculated as: Assuming that t l for l = 1, . . . , L are different from one another, with K > (U + 1)L and N M > (U + 1)L, it can be noted from the definition of matrix V that the dimension of its column space is equal to (U + 1)L with a full rank property. It follows that the dimension of the signal subspace is equal to the rank of R a . In the following we will show that R a is a full rank matrix of rank equal to (U + 1)L, which implies that the dimension of the signal subspace is (U + 1)L.
Proof Elements. Consider one path (L = 1) of mean delay t 1 ; for the sake of simplification, notations are abbreviated as and a (n,m) is a (U + 1)-length random column vector which is replaced by a new vector a defined as: under the following assumptions: . . , P} are independent real random variables having the same distribution. • α p and τ p are independent.
It follows that: For uniformly distributed τ p between [− s 2 , s 2 ], we have: According to the elements of matrix E[aa H ], it follows that its determinant is different from zero, (det(E[aa H ]) = 0) and as it is a square matrix, then it is a full rank matrix with rank equal to U + 1.

Channel Delay Estimation
The model in Equation (14) suggests the possibility of using subspace-based parameter estimation techniques. In the existing works with an exact common support [9,13], subspace methods such as multiple signal classification (MUSIC) [23] and estimation of signal parameter via rotational invariance techniques (ESPRIT) [24] are applied directly to estimate the delays. Although a non-strictly exact common support model is also taken into account in [9], the conventional ESPRIT and MUSIC methods are used with this model as well. Herein, in accordance with the delay spreading phenomenon, the channel model is reformulated, and a new subspace-based method is proposed to cope with such a model. The principle of subspace-based estimation methods is to find the parameters characterizing V that best match the signal subspace, which can be estimated from the measurements, based on the eigendecomposition of the observation covariance matrix. The problem of estimating the correct dimension of signal or orthogonal subspaces is of great importance in such methods, especially in scattering situations where this problem becomes more complicated.

Asymptotic Analysis of the Signal Subspace Dimension
In the exact common support case with no delay spreading, the noise-free covariance matrix R Y is of rank L. However, according to the Uth-order Taylor-approximated model (6), due to the delay spreading, the rank of R Y is increased to (U + 1)L, leading to a larger dimension of the signal subspace. More generally, considering the exact model (Equation (4)) with additive noise, most of the energy of the signal will be concentrated in few eigenvalues of the covariance matrix depending on the standard deviation of the delay spreading, and on the level of noise. Due to these factors, different dimensions of the signal subspace can be obtained. It follows that determining the appropriate value of U is of great importance, and it should be chosen carefully depending on these factors in order to estimate the so-called effective dimension of the signal subspace. Consequently, based on the definition of matrix V, an appropriate parametric expression of vectors can be chosen in the channel path delays estimation procedure.

Signal Subspace Tracking
From the generalized model in (6), the dimension of the signal subspace is (U + 1)L, where U is the order of the Taylor expansion that should be estimated to define the effective dimension of the signal subspace. This can be done through analyzing the eigenvalues of the estimated covariance matrix. The MDL and the Akaike information criterion (AIC) [25] are the two most famous information theoretic criteria for estimating the number of sources impinging on an array of sensors [26]. Such methods test the eigenvalues of the estimated covariance matrix where the number of sources is estimated based on minimizing a derived criterion. In general, the MDL is preferred as it is considered to have better performance than AIC [27]. The MDL provides better performance in the presence of spatially and temporally white noise [28]. The AIC provides better estimation in low signal-to-noise ratio (SNR) and small observation sample conditions, however this estimator is not asymptotically consistent. The MDL provides better estimation in the case of large numbers of observation samples, and it is asymptotically consistent. Therefore, the MDL criterion is chosen in this work to estimate the effective dimension of the signal subspace. In fact, the MDL criterion will provide (U + 2)L as the dimension of the signal subspace. Despite the fact that the error term R U (τ (n,m) l p ) in (5) is assumed to have no effect in describing the effective signal subspace, it will be detected by the MDL criterion as it is not equal to zero. Hence, during the signal subspace tracking process, the effective dimension of the signal subspace is estimated by subtracting L from the value provided by the MDL criterion. In fact, for L > 1, each of the L paths is assumed to occupy the same degrees of freedom (U + 1) in the signal subspace, which should give rise to an overall subspace dimension equal to (U + 1)L, and hence (U + 2)L would be provided by the MDL criterion. However in practice, in the presence of noise, the MDL algorithm will sometimes provide an integer in the set {(U + 1)L + 1, (U + 1)L + 2, . . . (U + 2)L}. The proposed solution for this case is to estimate U according to the following rule: where MDLV is the value obtained by the MDL criterion and . is the floor function.

Mean Path Delays Estimation
The covariance matrix of Y is estimated as in (15). As shown before, the delay spreading gives rise to an increase in the signal subspace dimension to (U + 1)L, where the value of U changes according to the SNR as well as the standard deviation of the delay spreading. It follows that the eigenvectors associated with the (U + 1)L largest eigenvalues define the signal subspace. Hence, the eigendecomposition of the estimated covariance matrix provides the noise subspace matrix U n , where U n is composed of the K − (U + 1)L eigenvectors associated with the K − (U + 1)L smallest eigenvalues of R Y .
For the case of no delay spreading, the conventional MUSIC can be applied to estimate the delays where the vector v(t) is projected on the noise subspace. In the case of delay spreading, as observed in Equation (6), vectors v (u) (t) for u = 0 . . . , U define the signal subspace. A new cost function is then proposed, such that vectors v (u) (t) are jointly projected onto the estimated noise subspace of dimension K − (U + 1)L. Hence, the following cost function can be defined: where The derived cost function can be written in the following form: where The proposed cost function turns to be the conventional MUSIC when setting U = 0, which is the case for very small delay spreading (where it tends to be neglected). The proposed technique can be seen as an extension of MUSIC in the delay spreading case.
The algorithm used for the estimation is summarized as follows:

Algorithm of the proposed method
Number of delays L is assumed to be known. 1. Collect the N M DFT-domain vectors to build matrix Y (14). 2. Estimate the covariance matrix R Y (15). 3. Apply the eigenvalue decomposition on R Y . 4. Estimate the effective dimension of the signal subspace using the MDL criterion combined with the proposed rule (22). 5. Construct matrix U n . 6. For each value of t, calculate the cost function P(t) (25). Then, the mean path delays t l are estimated by the positions of the L peaks of P(t).
The proposed cost function shows no major difference from the conventional MUSIC in terms of the computational complexity as they share the same operations such as the FFT, R Y matrix construction, eigenvalue decomposition, and the spectral searching. The proposed cost function requires a slight increase in the computational complexity due to the addition of the terms corresponding to the delay spreading, which is the projection of v (u) (t) for u = 0 . . . , U on the estimated noise subspace.

Simulation Results
Simulations are carried out for a 24 × 24 MIMO system but similar conclusions can be obtained with other MIMO system configurations. The proposed method can be used for either smaller or larger scale of MIMO systems. The DFT-domain data are generated according to Equation (4) with K = 64 Fourier coefficients. The delay deviations τ (n,m) l p of multi-rays within each path l for each (n, m) transmit-receive antenna pair are generated according to a uniform distribution and centered at t l .
The corresponding gains α (n,m) l p of the multi-rays are modeled as complex Gaussian random variables generated in a way that the effective power of each path is normalized, such that ∑ P l p=1 E[|α (n,m) l p | 2 ] = 1, ∀{l, n, m}. The added noise is modeled as complex white Gaussian noise. The estimation performance is assessed from Q = 500 independent simulations. The root mean square error (RMSE) of mean path delay estimation is calculated as , wheret l (i) is the estimated mean delay for the ith experiment and t l is the true mean delay. A unique delay group L = 1 is firstly considered. Figure 1 shows the comparison of RMSE of the mean delay estimation versus SNR of the proposed cost function for different values of U and MUSIC.  Figure 1 shows that the proposed cost function (assuming U = 1) provides a better performance of the mean path delay estimation than the conventional MUSIC.
According to the derivation of the model, in the noise-free case, it is better to choose higher values for U to benefit from a better approximation, where small details can be taken into account. On the other hand, with noise-contaminated data, the highest order terms in Taylor expansion (Equation (5)) may be too small with respect to the noise power. Therefore, due to the property of Taylor expansion, it can be noted from Equations (5) and (6), that: Hence, for U = 2 the elements in vector a (n,m) l,2 v (2) (t l ) will tend to be quite small. Thus, for a high level of noise, the norm of the vector a (n,m) l,2 v (2) (t l ) is very small compared to noise. Hence, this part of the signal is dominated by the noise, and it would be better not to consider this part of the signal to represent the signal subspace. For this reason, for low SNR, as can be observed in Figure 1, it would be better to consider U = 1 (K − 2L as noise subspace). However as SNR increases, it would be worth considering vector v (2) (τ l ) as a part of the signal subspace, and then U = 2 (K − 3L as noise subspace); this can be noted from the figure for SNR 10 dB. From Figure 2, we can observe that for very small delay spreading with relatively low SNR (SNR = 5 dB), it would be better to consider that the dimension of the signal subspace is one (U = 0), for the same reason discussed before, which is the situation of the conventional MUSIC. However, as delay spreading increases, the proposed cost function with U = 1 outperforms MUSIC (U = 0). It can be also noted that as the delay spreading increases, the proposed cost function for U = 3 would tend to outperform the proposed cost function for U = 2. The increase in the standard deviation of the delay spreading will give rise to an increase in a (n,m) l,u v (u) (t l ). Then, this vector will have greater impact in formulating the signal subspace, hence assuming higher value of U.
We can also observe that as the delay spreading increases, the estimation is less accurate, which is probably due to the approximation error of Taylor expansion.

Select U According to the MDL Criterion
To determine the most appropriate value of U, the MDL criterion is applied on the eigenvalues of R Y . The calculated mean value of the MDL criterion (mean MDLV) versus delay spreading is plotted in Figure 3 to show the performance of the MDL criterion.
However, it can be noted that the different flat regions are not always clearly distinguished. For example, when the delay spreading is within [−0.022, +0.022]T, from the observed mean MDLV values, it seems that the MDL criterion is providing four and five values in different realizations. This shows that the MDL criterion, influenced by the random nature of scattering and the relatively high noise level, is not always a foolproof indicator. In Figure 4, the RMSE of mean delay estimation of the proposed subspace tracking based method is shown. For each realization of received data, the decision about U is obtained from the rule in Equation (22), which is used to estimate the effective dimension of the signal subspace. Then, the proposed cost function (Equation (25)) is applied.  As shown in Figure 4, the proposed method allows for optimal selection of parameter U and provides the best mean path delay estimation. However some minor failures of the MDL criterion in systematically estimating the optimal effective dimension of the signal subspace can be noted. Figures 5 and 6 show the RMSE of the mean delay estimation and the mean MDLV versus delay spreading, respectively for SNR = 15 dB.  The obtained results show that the MDL criterion provides better estimations of the different effective dimensions of the signal subspace for different delay spreading intervals, and the different regions are well distinguished, leading to an improvement in the mean delay estimation performance of the proposed subspace tracking-based method. As mentioned before, the increase in the RMSE with respect to the delay spreading may be due to the approximation error in Taylor expansion. However, as shown in Figures 4 and 5, this increase is less significant for the proposed subspace tracking-based method; this is the main advantage of tracking the effective dimension of the signal subspace that changes according to the value of the standard deviation of the delay spreading and the noise level. Figure 7 shows the RMSE of mean delay estimation versus SNR of the proposed subspace tracking based method in comparison with MUSIC for different SNR values. The obtained results show a moderate improvement in the estimation performance as SNR increases. In fact as SNR increases, better estimation of covariance matrix is obtained, and hence there is better estimation of signal or noise subspaces.
As one path delay (L = 1) is considered before for the sake of simplification, Figures 8 and 9 show the RMSE of mean path delay estimation for L = 2 and L = 3 delays, respectively, for SNR = 15 dB.

Conclusions
A parametric subspace-based estimation method is proposed in this paper to deal with a non-strictly exact common support MIMO channel model that exhibits delay spreading due to scattering. A parametric model of the MIMO channel is first derived. Then, a subspace-based method is developed based on the tracking of the effective dimension of the signal subspace, which depends on the channel features. The proposed subspace tracking-based method is applied to estimate the channel mean delays. The proposed method shows better performance in comparison to the conventional methods, and allows estimating channel mean path delays with more accuracy, leading to a better MIMO channel path delay estimation in scattering environments.

Conflicts of Interest:
The authors declare no conflict of interest.

Abbreviations
The following abbreviations are used in this manuscript: