Nội dung

Hindawi Publishing Corporation EURASIP Journal on Advances in Signal Processing Volume 2007, Article ID 49393, 13 pages doi:10.1155/2007/49393 Research Article New Approaches for Channel Prediction Based on Sinusoidal Modeling Ming Chen,1 Torbjörn Ekman,2 and Mats Viberg1 1 Department of Signals and Systems, Chalmers University of Technology, SE 412 96 Göteborg, Sweden of Electronics and Telecommunications, Norwegian Institute of Science and Technology, NO-7491 Trondheim, Norway 2 Department Received 4 December 2005; Revised 4 April 2006; Accepted 30 April 2006 Recommended by Kostas Berberidis Long-range channel prediction is considered to be one of the most important enabling technologies to future wireless communication systems. The prediction of Rayleigh fading channels is studied in the frame of sinusoidal modeling in this paper. A stochastic sinusoidal model to represent a Rayleigh fading channel is proposed. Three different predictors based on the statistical sinusoidal model are proposed. These methods outperform the standard linear predictor (LP) in Monte Carlo simulations, but underperform with real measurement data, probably due to nonstationary model parameters. To mitigate these modeling errors, a joint moving average and sinusoidal (JMAS) prediction model and the associated joint least-squares (LS) predictor are proposed. It combines the sinusoidal model with an LP to handle unmodeled dynamics in the signal. The joint LS predictor outperforms all the other sinusoidal LMMSE predictors in suburban environments, but still performs slightly worse than the standard LP in urban environments. Copyright © 2007 Ming Chen et al. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited. 1. INTRODUCTION Link adaption techniques, such as multiuser diversity, adaptive modulation and coding, and fast scheduling hold great promise to improve spectrum efficiency. However, the improvement on the system capacity depends heavily on the predictability of the short-term channel fades [1, 2]. Extensive studies on this topic were made during the last several years by different researchers [1], [3–10]. It was found that a prediction horizon corresponding to a distance of half a wavelength traveled by the mobile is considered challenging [11]. In this paper, we assume the availability of a vector y = [y(t), y(t − 1), . . . , y(t − N + 1)]T containing the N channel observations (successive estimates of a particular channel coefficient) y = h + e, Figure 1, where the time index t = 0 and the length of the observation interval is N = 100. The published channel predictors are divided into two categories, which can be categorized as model-free predictors and model-based predictors, respectively. 1.1. Model-free channel predictors The first category is essentially the class of linear predictors (LP), where the channel coefficient is predicted as a weighted sum of the previous channel observations. A dth order LP of h(t + L), where d < N, is  + L) = h(t d −1 k=0 βk y(t − k) = βTd yd = flH y, (2) (1) where where h = [h(t), h(t − 1), . . . , h(t − N +1)]T is the time evolution of the true channel, and e = [e(t), e(t − 1), . . . , e(t − N + 1)]T is the additive estimation errors with Gaussian distribution CN (0, σe2 IN ), where IN is the N × N identity matrix. The value h(t + L) is to be predicted from observations y, where L is the prediction horizon. Such a scenario is presented in  flH = βTd 01xN −d  (3) is the coefficient vector of the LP. A large number of algorithms can be found in the literature to estimate β under various optimization criteria, for 2 EURASIP Journal on Advances in Signal Processing 25 where   T  Amplitude (dB) T φ = φ1 , . . . , φ p , ρ = ρ1 , . . . , ρ p , T ω = ω1 , . . . , ω p . 30 h(L) 35 These model parameters are assumed to be stationary over the observation interval. In these methods, the first step is to compute parameter estimates of ρ, φ, and ω. Then, the straightforward prediction of h(t + L) is  + L) = h(t y =h+e 40 100 80 (10) p   ρi e j φi e j ω i (t+L) . (11) i=1 60 40 20 Index of channel observations Let si = ρi e jφi be the complex amplitude. Then, (7) can be written as 0 Figure 1: Prediction of Rayleigh fading channel. y(t) = p  si e jωi t + e(t), (12) i=1 and the prediction of h(t + L) is instance, the least squares (LS) estimate, which is  β d = YH ls Yls where −1 † YH ls yls = Yls yls ,  yls = y(t), y(t − 1), . . . , y(t − N + L + d) (4) T  + L) = h(t (5) ⎤ · · · y(t − L − d+1) ··· y(t − L − d) ⎥ ⎥ ⎥. .. .. ⎥ . ⎦ . y(t − N +d) y(t − N +d − 1) · · · y(t − N +1) y(t − L) ⎢ y(t − L − 1) ⎢ Yls = ⎢ .. ⎢ ⎣ . y(t − L − 1) y(t − L − 2) .. . (6) 1.2. Model-based channel predictors In the second category, the channel over a short observation interval is modeled as superimposed deterministic complex sinusoids, which correspond to the Doppler frequencies as in the Jakes model [12], y(t) = h(t) + e(t) = p  ρi e jφi e jωi t + e(t), (7) i=1 where p is the number of sinusoids, ρi , φi , and ωi are the real amplitude, phase, and the Doppler frequency associated with the ith path, respectively, and e(t) is the additive complex Gaussian estimation error, e(t) ∼ CN (0, σe2 ). Assuming there are no moving reflecting objects, the ith Doppler frequency is ωi = 2πv cos θi , λ (8) where v is the velocity of the mobile terminal, and θi is the angle between the ith impinging ray and the mobile heading. In this model, there are (3p + 1) parameters to be estimated. They are  T ψ = ρT , φT , ωT , σe2 ,  si e j ω i (t+L) . (13) i=1 In vector form, the N observations can be represented as and the matrix Yls is a Hankel matrix, which is ⎡ p  (9) y = As + e, (14) where   A = a1 , . . . , a p ,  T ak = e jωk t , . . . , e jωk (t−N+1) ,  T s = s1 , . . . , s p . (15) (16) A thorough performance evaluation of linear predictions and channel predictions based on deterministic sinusoidal modeling were made by simulations and real world data in [6]. According to the reported results, the LP outperforms the deterministic sinusoidal modeling based methods. All these studies were performed in single-in single-out (SISO) systems. Later, similar results were reported in multiple-in multiple-out (MIMO) scenarios in [13]. It was also claimed that the calculation complexity of LP increases exponentially with the increase of the dimensions of the MIMO systems. This makes LP costly for high dimension MIMO channel prediction. This is one of the motivations to pursue channel prediction methods based on sinusoidal modeling, where a fewer number of model parameters is expected than the LP in MIMO scenarios. Beside these prediction methods, a nonlinear prediction of mobile radio channel using the MARS modeling was studied in [14]. In this study, the channel prediction is studied in the frame of sinusoidal modeling of the fading channel. First, a statistical sinusoidal model of the Rayleigh fading channel is proposed. Based on this model, three sinusoidal LMMSE channel predictors are given. Later a joint moving average and sinusoidal (JMAS) channel predictor is proposed, which alleviates the influence of the nonsinusoidal model components and leads to the Joint LS predictors. Ming Chen et al. 3 In this paper, the statistical sinusoidal model and the extensions into single-in multiple-out (SIMO) and MIMO systems are presented in Section 2. Three sinusoidal LMMSE predictors are given in Section 3. The JMAS prediction and Joint LS predictor are presented in Section 4. Section 5 presents an analysis of the measurement data and the stationarity properties of model parameters. The performance evaluations of these sinusoidal modeling-based methods by Monte Carlo simulations and real measurement data are presented in Section 6. Section 7 contains the conclusions. 2. STATISTICAL SINUSOIDAL MODELING In SISO scenarios, the stochastic sinusoidal modeling of radio channel has the same form as the deterministic sinusoidal model in (14), but the complex amplitudes of the sinusoids s are modeled as zero-mean random variables with covariance matrix, E[ssH ] = σs2 I p . This extension leads to a number of LMMSE channel predictors, which are presented in the following section. The statistical sinusoidal model for a SIMO system, with one transmit antenna and m receive antennas, is Y = AS + E, (17) where the matrices H are (Nm × n) and S is (pm × n), respectively. The power matrix is E[SS H ] = nσs2 I pm , where we have assumed, again, that the amplitudes are independent with equal mean power σs2 . Including channel estimation errors, the observed channel based on sinusoidal modeling for MIMO system becomes Y = AS + E . (22) Note that we have assumed here that all mn subchannels share the same Doppler frequencies. The extension to the general case is obvious, but introduces considerably more unknown parameters to be estimated. In the following sections, we restrict our study within SISO scenarios. 3. LMMSE PREDICTORS In this paper we are not concerned with the frequency estimation problem, but focus on prediction only. A number of successful frequency estimation methods have been published in the literature [16–19]. The class of subspace methods provides high-performance estimators at low cost, and of these we have chosen the unitary ESPRIT method due to [20]. where  T 3.1. Y = yT (t), yT (t − 1), . . . , yT (t − N + 1) ,  S = sT1 , sT2 , . . . , sTp  T (18) , T E = eT (t), eT (t − 1), . . . , eT (t − N + 1) , where the m-vector y(t) is the channel observations at the sensor array at time t. Matrix S is p × m, and the m-vector si is one realization of the complex random amplitude vector associated to the ith path. Assume these sinusoids have equal mean power σs2 , and E[SSH ] = mσs2 I p . The m-vector e(t) is the additive noise with covariance matrix E[eH (k)e(l)] = σe2 Im δk,l . Note that we have assumed that the angle of arrivals, θ, are identical at different antenna elements, but with independent amplitudes. In practice, inaccurate calibration of the array might give rise to such uncorrelated amplitudes. To derive the signal model for a MIMO system with n transmit antennas and m receive antennas, first let the SIMO channel, associated to the kth transmit antenna, be Hk = ASk , (19) where A and Sk are the same as in (17), except the subscripts. Vectorizing Hk , we have         hk = vec Hk = vec ASk = Im ⊗ A vec Sk ,   = Im ⊗ A sk ,   Assume that the frequency estimates ω  = [ω 1 , . . . , ω p ] are given, the LMMSE estimate of s is then [4]   H −1  H  A  + αI −1 A  H y = Rreg s = A A y, (23)  is defined as A, but using the estimated frequencies where A ω  k instead of the true frequencies, and α = σe2 /σs2 = 1/SNR. This can also be interpreted as a regularized LS estimator.  is Moreover, the LMMSE prediction of h(t + L) (given that ω assumed to be correct) is given by    + L) = e j ω 1 (t+L) , . . . e j ωp (t+L) s h(t −1  H = a(L)H Rreg A y = fcH y,  (24)  where a(L)H = e j ω 1 (t+L) , . . . , e j ω p (t+L) , −1  H fcH = a(L)H Rreg A (25) is the conditional prediction filter. If the complex amplitudes, s, are modeled as Gaussian random variables, CN (0, σs2 IN ), the statistical sinusoidal signal model (14) is termed Bayesian general linear model in [21], and the conditional LMMSE predictor is also the conditional MMSE predictor [10]. (20) 3.2. where ⊗ is the Kronecker product and vec(·) is the vectorization operation [15]. Stacking hk , the sinusoidal signal model for the MIMO channel is  Conditional LMMSE channel predictor  H = h1 , h2 , . . . , hn = Im ⊗ A][s1 , s2 , . . . , sn = AS, (21) Adjusted conditional LMMSE channel predictor Above, the frequency estimates ω  were regarded as exact. In  is estimated by some method, such as unitary ESpractice, ω PRIT [20], and subject to errors. The error variance, in theory, is determined by the SNR and the number of observations [21]. Such a problem was investigated in [22, 23]. Define the residual signal =[(t), . . . , (t − N + 1)]T as the part 4 EURASIP Journal on Advances in Signal Processing of the observations which cannot be accurately reconstructed by the estimated sinusoidal model, The expectation over the frequency estimation error can be expressed as   = y − h, E e j(ω k +Δωk )(n−m) = e j ω k (n−m) E e jΔωk (n−m)  (26)    = e j ω k (n−m) ΦΔωk (n − m), (36)  = A  s. In the imperfect modeling case, a low-order where h LP is required to predict the colored residuals (t + L), which contain less signal dynamics than y. So the adjusted conditional LMMSE predictor is where ΦΔωk (τ) is the characteristic function of the frequency estimation error Δωk . If we assume the frequency errors to be Gaussian distributed, then  + L) + (t + L). hadj (t + L) = h(t ΦΔωk (n − m) = e−(σΔωk /2)(n−m) . (27) The predictor (t + L) is computed based on past residuals as (t + L) = d −1 βk (t − k). 2 (t + L) =  βTd   (28) −1   reg 0TN −d IN − AR A  H y = fH y,   hadj (t + L) = fcH + fH y = faH y, (38) The (m, n)th element of the covariance matrix is thus obtained as  Rhh  mn = p  2 σs2 e j ω k (n−m) e−(σΔωk /2)(n−m) . k=1  Rhh   mn = σs2 p   2 e j ω k (n−m) e−(σΔω /2)(n−m) . [Γ]mn = e−(σΔω /2)(n−m) . Another approach to combat the frequency estimation error is motivated by introducing the frequency estimation uncertainty into the conditional LMMSE predictions. We can  k and variance model ωk as a random variable with mean ω 2 E[Δωk2 ] = σΔω , k  k + Δωk . ωk = ω (31) The classical form of the LMMSE predictor is given by the Wiener filter [21],  + L) = Rhy R−1 y, h(t (32) yy E    (33) R y y = E yy  = Rhh + Ree . Under the assumption that the amplitudes and frequency estimation errors are independent stochastic variables, we can write the (m, n)th element of Rhh as  Rhh  mn =E  AssH AH   k=1 (42) where  denotes the Hadamard product. Similarly, the cross covariance Rhy is given by Rhy p     2  H H H =E  a(L) ss A = E sk   k=1 ·E e j(ω  k +Δωk )L ,...,e j(ω  k +Δωk )(L+N −1)  (43) . Define the damping vector 2 2 2 2  (44) The cross covariance is then obtained as p   e j ω k L , . . . , e j ω k (L+N −1)  γ k=1 (45)  H  γ. = σs2  a(L)H A The unconditional LMMSE predictor is finally given by (just some variation)  + L) = Rhy R−1 y h(t yy mn p      = E |sk |2 E e j(ω k +Δωk )(n−m) . (41) A  H  Γ, Rhh = σs2 A Rhy = σs2 (34) (40) The covariance matrix can then be expressed as  The covariance matrix for N observations is H 2 γ = e(−σΔω /2)L , . . . , e(−σΔω /2)(L+N −1) .   h(t + L)  σh2 Rhy h(t + L)H yH = . y R yh R y y l  2 k=1 2  (39) Define the (N × N) damping matrix Γ by 3.3. Unconditional LMMSE predictor where 2 2 For simplicity, let all the frequency errors be IID with σΔω = k 2 σΔω . We then have (30) where faH is the adjusted conditional LMMSE prediction filter.  E |sk |2 = σs2 . (29) where 0N −d is an (N − d) zero vector, and fH is the residual predictor. Then, the adjusted conditional LMMSE predictor (27) can be expressed as (37) The expectation over the ensemble of amplitudes is just the variance k=0 The coefficients βd = [β0 , . . . , βd−1 ]T are computed by solving an LS problem. To explicitly show the dependency of (t + L) on the channel observation y, we rewrite (28) as 2 (35)  H   A H  γ A   Γ + αIN −1 y =  a(L)H A = fuH y, (46) Ming Chen et al. 5 where   A H  γ A  H  Γ + αIN fuH = a(L)H A To take into account the prediction model error, (50) can be rewritten as −1 (47) is the unconditional prediction filter. Clearly, the influence of old observations is reduced by the damping matrix Γ and the damping vector γ, which are also dependent on the prediction range. The further ahead we are looking, the smaller the gain of the predictor is. The way the frequency error is taken into account can be interpreted as a convolution of the line spectrum of the signal with the error distribution. The filter design is thus done for distributed sources. The errors in the frequency estimates force the predictor to misstrust older data, as even a small frequency error can cause a large phase error if one waits long enough. As the special 2 2 case, when σΔω = 0, (46) degenerates to (24). When σΔω > 0, (46), although being LMMSE (ignoring the influence of sk on the estimate of ωk ), is not the MMSE due to the nonlinear dependency on the frequency estimates. 4.  + L) = y T β +  h(t a(L)H s + e(t + L), d d (48) where ω  in  a(L) contains only the frequencies of the stationary sinusoidal signals, which are determined by some model selection method. With this model, the channel is divided into two parts. One part contains the consistent sinusoidal signals, such as the direct rays in LOS scenarios. The second part captures all the remaining power in the observations. This prediction model is interpreted as the joint moving average and sinusoidal (JMAS) predictor. Note that the name refers to the predictor, but not to the signal model. Selecting which sinusoids to use in the predictor is not a trivial problem, and deserves further study. A possibility is to use frequency tracking and apply some stationarity measure for each frequency track. Due to the linearity of the autoregressive and sinusoidal bases, the model parameters, θ = [βTd sT ]T , can be estimated jointly. Given ω,  (48) can be rewritten as    β d  s . (49) In vector form, we have    N −L−d+1 = Yls A  J β = Hθ, h (50) where T   − 1), . . . , h(t  + L − N + d) ,   N −L−d+1 = h(t), h h(t    J =  A aJ1 , aJ2 , . . . , aJ p ,  aJk = e j ω k (t+L) , e j ω k (t+L−1) , . . . , e j ω k   (t+2L−N+d) T (51) . −1 θ = HH H HH yls . (53) The channel prediction is   H   + L) = y T  h(t θ, d a(L) (54) which is named Joint LS predictor. In the terminology of the  is previous sections, it is a conditional LS predictor, since ω assumed to be given and s is modeled as deterministic. Finally, substitute (53) into (54) to obtain     a(L)H = ydT  In order to represent signals with more general spectral characteristics, let the estimate of h(t + L) be modeled by (52) where u = [u(t), u(t − 1), . . . , u(t − N +L+d)]T with Gaussian distribution CN (0, σu2 IN −L−d+1 ). Then, the LS estimate of θ is H  + L) = y T  h(t d a(L) JMAS PREDICTION MODEL AND JOINT LS PREDICTOR H  + L) = y T  h(t d a(L)  N −L−d+1 + u yls = h = Hθ + u, −1 HH H HH yN −L−d+1  (HH H)−1 HH 0M y = fJH y, (55) where the 0M is a zero matrix of conformable dimensions, and  fJH = ydT a(L)H   −1 HH H HH 0M  (56) is the joint LS prediction filter. Note that this is a nonlinear predictor, since H depends on y. The difference between the joint LS method and the adjusted conditional LMMSE predictor in Section 3.2 can be seen by rewriting (26) as  H −1 H  A  A   y  =y−A A   H −1 H   A  A   y = Π⊥ = IN − A A A y. (57) So that (27) becomes     β T H d  hadj (t + L) = d a(L) =   s 1Td T 0TN −d ⊥  ΠA y T   a(L)H (58) θ, where 1Td is a d-vector with all ones. It can be seen that the autoregressive bases are orthogonal to the sinusoidal bases in (58), but they are not in (49). In general, the total order of the joint LS prediction, r = p + d, can be determined by some classical criterion, that is, AIC [24] and MDL [25]. So a wide sense definition of model selection includes the selection of the stationary sinusoidal part of the signals, p, and the order of the LP, d, or both at once. In this paper, the model selection refers to the selection of the strong and stationary sinusoidal signals only. In practice, it is probably a wise idea to overestimate the rank, and estimate unnecessarily many sinusoids. But only a subset of these is used in the predictor. The tricky part is to know which ones to keep. With only one data set, one would probably go for the ones with the largest estimated 6 EURASIP Journal on Advances in Signal Processing ANALYSIS OF MEASUREMENT DATA 5.1. Measurement data A measurement campaign was performed in the urban area of Stockholm, and the suburban area, Kista, a few kilometers north of central Stockholm. The measurements were performed at 1880 MHz, and the bandwidth is 5 MHz. In total, 45 and 31 effective measurements were performed at different urban and suburban locations, respectively. The velocities of the mobile station were between 30 to 90 km per hour during the measurements, except for a few stand still cases. Each measurement records channel sounding data over 156.4 ms and contains 1430 repetitions of channel impulse responses. This gives rise to the channel update rate of 9.1 KHz. Each channel impulse response (or power delay profile) is described by 120 taps with a sampling frequency of 6.4 MHz. The time interval between two neighboring taps is 0.156 μs. In this study, the strongest channel tap in each measurement is used in the performance evaluations of the proposed channel prediction methods, which can be considered as a narrow band Rayleigh fading channel. Obviously, the performance of the sinusoidal modelingbased channel prediction methods depends heavily on the stationarity properties of model parameters in the observation window and the prediction horizon. In this section, the stationarity of the model parameters is investigated by means of the short-term MUSIC pseudospectrum (MPS) and model parameter tracking. Frequency (rad/s) 5. 0.1 5 4 0.05 3 0 2 1 0.05 0 1 0.1 10 20 30 Number of blocks 40 Figure 2: Example of a MUSIC pseudospectrum of a measurement in an urban area. A number of high-power frequency bins can be observed. They are distributed on both the positive and negative sides of the spectrum. But the relative power between these bins and the remaining frequency bins is low. Those bins close to the positive and negative maximum Doppler frequencies are more consistent over the measurements than the others. These bins have blurry edges. 0.1 4 Frequency (rad/s) amplitude. With several data sets, one could try to evaluate the stationarity of the various frequency estimates, but this requires a clever sorting (data association) to create frequency tracks. 0.05 3 2 0 1 0.05 0 5.2. MUSIC pseudospectrum 0.1 The MUSIC pseudospectrum (MPS) is defined as [16] pmu (ω) = a(ω)H a(ω) , a(ω)H Π⊥ a(ω)   ω ∈ − fmax , fmax , Un UH n,  R y y = Us Un = Us Λ s UH s   Λ s  Λn UH s UH n 20 30 Number of blocks 40 (59) where Π = and the columns of Un are the eigenvectors spanning the noise space of the covariance matrix of the data sequence y, R y y = E[yyH ]. It can be obtained by taking eigenvalue decomposition of R y y as ⊥ 1 10  Figure 3: Example of a MUSIC pseudospectrum of a measurement in a suburban area. A fewer number of high-power frequency bins are observed compared to Figure 2. The locations of these frequency bins are biased to the negative side in the spectrum. The edges of these high-power frequency bins are much sharper than those in Figure 2. The power in the strongest frequency bin is much higher than the others and is consistent over the whole measurement interval. (60) + Un Λ n U H n, where a rank estimation algorithm is necessary, and the eigenvalues are ordered in nonincreasing order. In case of high SNR, MPS is not sensitive to under-estimated noise subspace dimension, while an overestimated noise space will attenuate certain weak frequency components. The a(ω) is the DFT vector associated with frequency ω, and fmax is the maximum Doppler frequency. The MPS is calculated over each segmented data block, where a sliding window is applied for data extraction. Examples of typical time varying MPSs in urban and suburban channels are presented in Figures 2 and 3, respectively. In these two measurements, the mobile speeds were around 30 km/h, the original data is downsampled by a factor of 10, and the window length (after downsampling) is 100, Ming Chen et al. which corresponds to a spatial observation interval of 6 λ. The numbers of sinusoids are set to be 8 and 7, respectively, which is assumed to be slightly overestimated. In these figures, the grayscale of the contours indicates the powers in the frequency bins. In Figure 2, the MPS from an urban measurement is given, where a number of frequency bins with deeper grayscale appear in the spectrum. A number of relatively high-power and stationary frequency bins, close to the positive and negative maximum Doppler frequencies, can be observed. Their edges are blurry, which implies that these bins are contributed by a cluster of (or more than one) close, but separated reflectors. This is coincident with our simulation results. The remaining frequency bins, close to zero frequency, are weaker and less stationary. This is due to that the low Doppler frequencies are associated with the impinging rays perpendicular to the mobile velocity. The observed frequencies and amplitudes belonging to these rays are more likely to experience large time variation over the observation interval compared to the maximum Doppler frequency components close to the direction or the inverse direction of the mobile velocity. Such a spectrum agrees with the typical Doppler spectrum in a rich scattering urban environment. These nonstationary frequency bins might render sinusoidal modeling-based channel prediction difficult. The MPS in Figure 3 is obtained from a measurement in a suburban area. It has much fewer high-power frequency bins compared to what is observed in urban measurement as given in Figure 2. The distribution of these frequency bins in the spectrum is biased to the negative side in this example. This is probably due to that the mobile station was moving away from the base station or the reflection objects during the measurement. The edges of these high-power frequency bins are much sharper than those in Figure 2. This indicates that each bin might contain only one pure sinusoid. This is also coincident with a typical Doppler spectrum in a suburban or rural area in LOS. The power in the strongest frequency bin is much higher than the others, and is consistent over the whole measurement interval. 7 measurement, as given in subplot (a) in Figure 5. The channel phase is linear in just a part of the measurement, which indicates that there might be just one dominant sinusoid in the linear part of the channel, but more comparable dominant sinusoids in the remaining part. In subplot (d), one amplitude is much higher than others, which is marked by an asterisk. Its frequency is close to −0.1 as seen in subplot (c). The amplitude of this ray is consistent over the observation interval. Besides this sinusoid, a number of sinusoids with well-separated frequencies are found in subplot (c), but their amplitudes are time varying and comparable. The variations of these amplitudes are much larger than that of the first frequency. The presence of these sinusoids might explain the large variation of the channel amplitudes and nonlinear phase in this measurement. In Figure 5, the parameter tracking using the suburban example measurement is given. Both the average and the variation of the channel amplitude are smaller than in the urban measurement, as seen in subplot (a) in Figure 4. The phase of the channel is quite linear over the whole measurement interval. In subplot (d) there are two sinusoids marked by a triangle and a circle having the strongest amplitudes among all the estimated sinusoids. The frequencies associated with these two sinusoids are extremely close, as seen in subplot (c). This may be due to over estimated number of sinusoids and indicates that the channel could be well modeled by just one dominant sinusoidal component. 6. 6.1. PERFORMANCE EVALUATIONS Measures of prediction performance The performance of the predictors is evaluated by Monte Carlo simulations and measurements. The normalized mean square error (NMSE) is adopted in this paper for performance evaluations. First, we define the normalized square error (NSE) in a single realization as  2 eNSE (t 2  + L) N h(t + L) − h(t + L) = . H h h (61) 5.3. Model parameter tracking Similar to the study of the MPSs, estimates of the model parameters (frequencies and amplitudes) are computed for each data block using the unitary ESPRIT and LS methods. The estimates of the frequencies and amplitudes using the same urban and suburban example measurement are plotted as a function of the number of data blocks in Figures 4 and 5, respectively. The numbers of sinusoids are assumed to be 8 and 7 for urban and suburban example measurements, respectively, as before. In these figures, the channel amplitudes and phases are given in subplots (a) and (b). The estimated frequencies and the associated amplitudes are given in subplots (c) and (d), where different markers are used for frequency and amplitude estimates associated with each sinusoid. In Figure 4, which is the urban case, the dynamic range of the channel amplitude is larger than the one in the suburban 2 (t + L) taken over Then, the NMSE is the mean value of eNSE the Monte Carlo simulations. However, an adjusted NMSE (ANMSE) is also introduced to get rid of the influence of a small amount of outliers which might ruin the whole averaged performance. So the outliers are dropped when the NSE of their power prediction is larger than 0.04, or in other words, when the prediction error of the complex amplitude is larger than 20% of the root mean square (RMS) channel amplitudes. The concept of level of confidence is also used in the performance evaluation. For example, when we say the NMSE of a channel predictor is −10 dB with the level of confidence of 95%, we mean that the average of the NSEs of the best 95% of the channel predictions from simulations or measurements is less than −10 dB, while the worst 5% cases are excluded. It is also worth noting that when the measurement data is used, the channel observation y(t + L) is used as h(t + L), 8 EURASIP Journal on Advances in Signal Processing 10 4 3 2 1 0 2 Channel phase (rad) Channel amplitude 3 0 2 4 0 50 100 Number of samples 0 50 100 Number of samples (b) 10 2 10 4 10 6 0.2 Estimated amplitudes Estimated frequencies (rad/s) (a) 0.1 0 0.1 0.2 0 10 20 30 Number of blocks 40 0 10 20 30 Number of blocks (c) 40 (d) Figure 4: Example of the time variation of model parameters in an urban area: (a) the channel amplitude, whose dynamic range is large in this example; (b) the phase of the channel, which is linear in a part of this measurement; (c) the estimated frequencies of the sinusoids in each block, the number of sinusoids is assumed to be 8; (d) the amplitudes of the sinusoids associated to the frequencies in (c). and the prediction error is normalized by the mean power of y, instead of h as given in (61). 6.2. Model selection based on SVD In this study, an ad hoc singular-value-decomposition(SVD-) based model selection method is proposed. First, the number of the consistent sinusoidal signals is estimated using the ith largest gradient method. It chooses the index of the descending ordered singular values which gives the ith largest gradient in linear scale. For example, an m-vector  σ = σ1 , . . . , σr −1 , σr , . . . , σm T (62) contains the singular values from the SVD of Yls , where L = 0 and d = m. The value of m should be larger than the number of sinusoids. These singular values are arranged in descending order. If (σr −1 − σr ) gives the ith largest gradient between adjacent singular values, σr is called the break point, and the index r is then selected as the number of consistent sinusoidal signals. In the second step, p frequencies are estimated from the data block, where p > r. The frequencies that have the r largest estimated amplitudes are, then, assumed to be consistent. The selection of i in this simple method is environment dependent. In a suburban area, most rays can be expected to be consistent. We could therefore choose a larger i. But in an urban area, this method might not fit, since some sinusoids are strong, but not stationary. This means that a small i should be selected for urban measurements to be conservative. In the following performance evaluations, we choose i = 1 and i = 3 for urban measurements and suburban measurements, respectively, which provide the best performance for the methods in Section 2. 6.3. Simulation setup The simulation parameter setting is as follows: (i) SISO scenario; (ii) spatial channel sampling interval = 0.1 λ (Δl = 0.1); (iii) prediction horizon = 0.5 λ (L = 5); Ming Chen et al. 10 9 4 3 3 2 Channel phase (rad) Channel amplitude 2.5 2 1.5 1 0 2 0.5 4 0 50 100 Number of samples 0 50 100 Number of samples (b) (a) Estimated amplitudes Estimated frequencies (rad/s) 0.2 0.1 0 0.1 0.2 0 10 20 30 Number of blocks 10 40 5 0 (c) 20 Number of blocks 40 (d) Figure 5: Example of the time variation of model parameters in a suburban area: (a) the channel amplitude, whose dynamic range is smaller than in Figure 4; (b) the phase of the channel, which is linear in the whole measurement interval; (c) the estimated frequencies of the sinusoids in each block, the number of sinusoids is assumed to be 7; (d) the amplitudes of the sinusoids associated to the frequencies in (c). (iv) number of sinusoids = 8 (p = 8, which is assumed to be known); (v) number of samples = 100 (N = 100); (vi) number of Monte Carlo simulations = 1000; (vii) order of LP is two times the number of sinusoids; (viii) in the combined method, the model selection is based on the powers of sinusoids. Specifically the sinusoids with |sk | > 2σe are predicted by the conditional LMMSE method, the residual is predicted using an LP with order 2 (d = 2); (ix) the frequency error variance is used as a design param2 = 1/N 3 . This is somewhat eter, and is taken to be σΔω arbitrary, but is motivated by the CRLB [21]; (x) the unitary ESPRIT method is used for frequency estimation. In Figure 6, the cumulative density function (CDF) of the NMSE is presented. All sinusoidal model-based predictors have light tails, but the LP has a heavy tail. The ANMSE of different predictors is given in Figure 7. It can be seen that the conditional LMMSE predictor with known frequency and the LP have the best and the worst performances, respectively. The conditional LMMSE predictor with estimated frequency has much better performance than the LP in the investigated scenarios. The unconditional LMMSE predictor performs slightly better than the conditional LMMSE predictor. The performance of the adjusted conditional LMMSE predictor approaches that of the conditional LMMSE predictor with known frequency with the increase of SNR. Note that, in these simulations, for a given carrier frequency, f , and velocity, v, the spatial sampling interval can be easily converted into time sampling interval as Δt = (Δl · c)/( f · v), where c is the speed of radio propagation. 6.4. Performance evaluation of LMMSE predictors using measured data The original data is downsampled by a downsampling ratio (DSR) of 10 to reduce the calculation load and increase the prediction horizon, which gives rise to a sampling frequency 10 EURASIP Journal on Advances in Signal Processing 0 1 5 10 NMSE (dB) Probability 0.8 0.6 0.4 15 20 25 30 0.2 35 0 40 0 10 20 Average SNR (dB) LMMSE, known ω Conditional LMMSE Unconditional LMMSE 30 40 ANMSE 1 10 2 10 3 10 4 10 5 10 6 0.4 0.6 Level of confidence LP Conditional LMMSE LP Adjusted Figure 6: Probability of channel prediction error less than 20% (L = 0.5 λ, N = 100, 8 rays). 10 0.2 0.8 1 Unconditional LMMSE Adjusted Figure 8: Performance evaluation of LMMSE prediction methods in an urban area. (iii) the SNRs used in the conditional LMMSE prediction methods are estimated from the associated power delay profiles; (iv) the order of the LP in the adjusted conditional LMMSE predictors is 2, which is fixed; (v) the order of the standard LP is set to be the same as the signal order estimated for the data block using model selection method; (vi) the 1st and 3rd largest gradients are used for the model selection based on SVD for urban and suburban environments, respectively; (vii) other settings are the same as in Section 6.3. 0 10 20 Average SNR (dB) LMMSE, known ω Conditional LMMSE Unconditional LMMSE 30 40 LP Adjusted Figure 7: Adjusted NMSE of channel prediction (L = 0.5 λ, N = 100, 8 rays). of 910 Hz. After the downsampling, the evaluation parameter settings are (i) number of samples = 100 (N = 100), which is corresponding to a measurement interval of 0.92 m (6 λ) or 2.75 m (18 λ) in distance at the speed of 30 km/h or 90 km/h, respectively; (ii) the prediction horizon L is 5, which is λ/3 or 1 λ at the above speeds and the frequency band of 1880 MHz; in time scale, it is a prediction of 5.5 ms into the future. With these settings there are 39 blocks of data in each measurement. In total, 1755 and 1209 blocks of data are obtained in urban and suburban areas, respectively. In Figure 8 the overall NMSEs of different prediction methods using all measurement data in an urban area are presented versus the level of confidence. It is meaningless to discuss the achievable NMSEs in this case due to the limited number of measurements. We put our attention on the relative performance instead. In this figure, the LP has the best performance. The conditional LMMSE predictor has the worst performance. The unconditional LMMSE predictor outperforms the conditional LMMSE predictor. The adjusted conditional LMMSE predictor outperforms the conditional LMMSE predictor and the unconditional LMMSE predictor, and is slightly worse than LP. In Figure 9, the results using suburban measurements are given. Similar results are observed. In these results, the performances of the conditional LMMSE and the unconditional LMMSE predictor are even closer compared to those in an urban area. This is mainly due to that the frequency separation is larger in suburban than in urban environments. The adjusted conditional LMMSE predictor outperforms the