ω(q): fω (x) = ln(10) 10(x+ΥMP, dB )/10 fΥMP 10(x+ΥMP, dB )/10 . 10 Some experimental measurements have also suggested Gaussian to be a good enough yet simple fit for the distribution of ω(q) [17]. We will take advantage of this Gaussian simplification later in our framework. As for the shadowing variable, log-normal is shown to be a good match for the distribution of ΥSH (q). Then, we have the following zero-mean Gaussian pdf 2 1 for the distribution of ν(q): fν (x) = √2πα e−x /2α , where α is the variance of the shadowing variations around path loss. Characterizing the spatial correlation of ω(q) and ν(q) is also considerably important for our model-based channel prediction framework. However, we do not attempt to predict the multipath component, ω(q), due to the fact that it typically decorrelates fast and that the form of its correlation function can change considerably, depending on the angle of arrival and position of the scatterers. Therefore, in our proposed framework we only predict the path loss and shadowing components of the channel. The impact of multipath will then appear in the characterization of the prediction error variance, as we shall see. As for the spatial correlation of shadowing, [18] characterizes an exponentially-decaying spatial correlation function, which is widely used: E ν(q1 )ν(q2 ) = α e− q1 −q2 /β , for q1 , q2 ∈ K where α denotes the shadowing power and the correlation distance, β, controls the spatial correlation of the channel [18]. For some examples of a time-varying environment, Oestges et al. model the slow temporal-variation of the channel as a zero-mean Gaussian variable with an exponential temporal correlation in the dB domain [27]. Thus, Eq. 1 can be extended to such time-varying cases by adding this additional variable to Eq. 1. In general, however, finding one model for characterizing the time-variations of different features in the environment is a challenging task and a subject of further studies. Such temporal variations can also be treated as disturbance in the prediction process. Next, we describe our proposed model-based channel prediction framework. Consider the case where a wireless channel to a fixed transmitter is sparsely sampled at positions Q = {q1 , q2 , · · · , qk } ⊂ K, in a given environment. These channel measurements can be gathered by one or a number of cooperative homogenous robots, equipped with identical receivers, making measurements along their trajectories. Let a region or an environment refer to an area over which the underlying channel parameters, such as α and β, can be considered constant. The four marked areas of Fig. 10 are examples of such regions. First, consider the case that all the k measurements belong to one region and that we are predicting the channel in the same region. We show how to relax this assumption later in this section. Let DQ and YQ = [y1 , · · · , yk ]T ∈ Rk denote the corresponding distance vector to the transmitter in dB and the vector of all the available channel measurements (in dB) respectively: DQ = 10 log10 ( q1 − qb ), · · · , 10 log10 ( qk − qb ) YQ = [y1 , · · · , yk ]T ∈ Rk . We have, YQ = 1k −DQ θ + ϑQ + ΩQ ,
HQ T

[KdB nPL ]T is the vector of the path loss parameters, ϑQ = T T ν1 , · · · , νk with νi = ν(qi ) and ΩQ = ω1 , · · · , ωk with ωi = ω(qi ), for i = 1, · · · , k. Based on the lognormal model for shadowing, ϑQ is a zero-mean Gaussian random vector with the covariance matrix RQ ∈ Rk×k , where RQ i,j = α e− qi −qj /β , for qi , qj ∈ Q. The term ΩQ denotes the impact of multipath fading in dB domain. As mentioned earlier, some empirical data have shown Gaussian to be a good match for the distribution of wi [17]. For instance, Fig. 2 compares the match of both Nakagami and lognormal to the distribution of multipath fading (ΥMP ) for a stationary section of our collected data of Fig. 1. As can be seen, Nakagami provides a considerably good match while lognormal can be acceptable, depending on the required accuracy. Thus, in order to facilitate the mathematical derivations in our prediction framework, we take wi to have a Gaussian distribution. In addition, multipath fading typically decorrelates considerably fast, making learning of its correlation function, based on sparse possibly non-localized samples, considerably challenging if not infeasible. There is also no one general function that can properly model its correlation in all the environments as its form depends heavily on the angle of arrival and position of the scatterers. While approaches based on the estimation of the power spectrum and linear prediction have been utilized to predict the immediate values of multipath, based on past observations, such approaches require dense sampling in order to capture correlated multipath samples. Finally, even if its correlation function is learned, it typically can not be taken advantage of, in the prediction framework, unless the location of the channel to be predicted is very close to the position of one of the available measurements. Thus, we take ΩQ to be an uncorrelated zero-mean Gaussian vector with the covariance of E ΩQ ΩT = σ 2 Ik×k , where Ik×k is a k × k identity Q matrix and σ 2 = E ω 2 (q) = 100
∞ 2 ∞ 0

log2 (x)fΥMP (x)dx − 10

and (2)

100 0 log10 (x)fΥMP (x)dx is the power of multipath fading (in dB domain). In other words, our framework does not attempt to predict the multipath component and assumes the worst case of uncorrelated multipath (worst from a prediction standpoint). The estimated variance of multipath then appears in our assessment of channel prediction error variance, as we shall see. Note, however, that this is only for the purpose of our modeling. When we show the performance of this framework, we use real measurements where the multipath component will have its natural distribution and correlation function. We then define ΞQ ϑQ + ΩQ , which is a zero-mean Gaussian vector with the covariance matrix of Rtot,Q RQ + σ 2 Ik×k . In our model-based probabilistic framework, we first need to estimate the parameters of the model (θ,α, β and σ 2 ) and then use these parameters to estimate the channel. Let fYQ (YQ |θ, α, β, σ 2 ) denote the conditional pdf of YQ , given the parameters θ, α, β and σ 2 . Under the assumption of independent multipath fading variables, Eq. 2 will result in the following: fYQ (YQ |θ, α, β, σ 2 ) = e− 2
1 T

YQ −HQ θ

αRnorm,Q (β)+σ2 Ik×k

−1

YQ −HQ θ 1/2

(2π)k/2 det αRnorm,Q (β) + σ 2 Ik×k

where 1k denotes the vector of ones with the length of k, θ =

MALMIRCHEGINI and MOSTOFI: ON THE SPATIAL PREDICTABILITY OF COMMUNICATION CHANNELS

967

8000

Probability Density Function

7000 6000 5000 4000 3000 2000 1000 0 −1000 0 1 2 3 4 5 6 −4 x 10

experimental data power distribution of Nakagami fading with m= 1.20 Lognormal fit

ML estimation of the channel parameters for the general case, where σ 2 = 0, is computationally complex. Therefore, we next devise a suboptimum but simpler estimation strategy. Let χ = α + σ 2 denote the sum of the shadowing and multipath powers. A Least Square (LS) estimation of θ and χ then results in:
−1 T T ˆ θLS = HQ HQ HQ YQ , (5) 2 1 T −1 T T χLS|θ=θLS = YQ Ik×k − HQ HQ HQ ˆ HQ YQ ˆ k 1 T −1 T T = YQ Ik×k − HQ HQ HQ HQ YQ , (6) k where HQ is full rank, except for the case where the samples are equally-distanced from the transmitter. Since such a special case is very low probable, we assume that HQ is full rank throughout the paper unless otherwise is stated. We refer to this suboptimal approach as LS throughout the paper. We next discuss a more practical but suboptimum strategy to estimate β. Let Il = (i, j) qi , qj ∈

Received power (mW)

Fig. 2: Comparison of Nakagami and lognormal for the distribution of small-scale fading.
1 where Rnorm,Q = α RQ denotes the normalized version of RQ . Next, we characterize the Maximum Likelihood (ML) estimation of the underlying channel parameters.

ˆ ˆ ˆ ˆ2 [θML , αML , βML , σML ] = argmaxθ,α,β,σ2 ln fYQ (YQ |θ, α, β, σ 2 ) = argminθ,α,β,σ2 YQ − HQ θ
T

αRnorm,Q (β) + σ 2 Ik×k ,
−1

−1

× YQ − HQ θ + ln det αRnorm,Q (β) + σ 2 Ik×k which results in:
T ˆ ˆ θML = HQ αML Rnorm,Q (βML ) + σML ˆ ˆ2 T ˆ × HQ αML Rnorm,Q (βML ) + σML ˆ2 ˆ −1

HQ

−1

Q such that qi − qj = l denote the pairs of points in Q which are located at distance l from each other. Let −1 T T HQ YQ represent the YQ,cent,LS = Ik×k − HQ HQ HQ centered version of the measurement vector, when path loss parameters are estimated using the LS estimator of Eq. 5. 1 Define rQ (l) ˆ (i,j)∈Il [YQ,cent,LS ]i [YQ,cent,LS ]j to be |Il | the numerical estimate of the spatial correlation function at distance l, where |.| represents the cardinality of the argument set and [.]i denotes the ith element of the argument vector. We have [ˆ LS , βLS ] = arg minα,β l∈LQ w(l) ln αe−l/β − α ˆ ln rQ (l) , where LQ = {l|0 < rQ (l) < χLS|θ=θLS } and ˆ ˆ ˆ ˆ w(l) can be chosen based on our assessment of the accuracy of the estimation of rQ (l). For instance, if we have very few pairs ˆ of measurements at a specific distance, then the weight should be smaller. Let LQ = {l1 , l2 , · · · , l|LQ | } denote an ordered set of all the possible distances among the measurement points. We have the following Least Square estimator of α and β: ln(ˆ LS ) α T T = (MLQ WLQ MLQ )−1 MLQ WLQ b where MLQ = 1 ˆLS ⎤ ⎤ ⎡ β ⎡ 1 −l1 ln rQ (l1 ) ˆ ⎥ ⎢. . ⎥,b = ⎢ . . ⎦ . = ⎦ and WLQ ⎣. ⎣ . . . k 2

YQ .

(3)

ˆ Finding a closed-form expression for αML , βML and σML , howˆ ˆ2 ever, is challenging. For the special case where Ωk is negligible, the ML estimation of channel parameters can be simplified to: ˆ θML,σ2 =0 =
T −1 ˆ HQ Rnorm,Q (βML,σ2 =0 )HQ −1 T −1 ˆ HQ Rnorm,Q (βML,σ2 =0 )YQ , −1 ˆ Rnorm,Q (βML,σ2 =0 )

αML,σ2 =0 = ˆ

1 ˆ YQ − HQ θML,σ2 =0 k ˆ × YQ − HQ θML,σ2 =0 , β T

ˆ βML,σ2 =0 = arg min

−1 T T YQ PQ,ML (β)Rnorm,Q (β)PQ,ML (β)YQ

× det Rnorm,Q (β) , where PQ,ML (β)
T −1 HQ HQ Rnorm,Q (β)HQ −1

(4) = Ik×k −

Under the ˆML,σ2 =0 is assumption that β is known, it can be shown that θ an unbiased estimator and achieves the Cramer-Rao bound. Furthermore, for large number of sampling points k, we can show that αML,σ2 =0 is unbiased and achieves the Cramer-Rao ˆ bound as well. We skipped the details of the proofs due to the space limitation. The ML estimator will therefore be our benchmark in the estimation of the channel parameters. As can be seen, in order to estimate θ and α, we first need to estimate β, which is challenging. Furthermore, finding the

T −1 HQ Rnorm,Q (β).

1 −l|LQ | ln rQ (l|LQ | ) ˆ ˆ2 ˆ diag[w(l1 ), · · · , w(l|LQ | )]. We then have, σLS = χLS|θ=θLS − ˆ αLS for the estimation of the multipath power (in dB domain). ˆ Note that the estimated values of the shadowing parameters ˆ ˆ should satisfy: 0 < αLS ≤ χLS|θ=θLS and βLS > 0. If due to ˆ ˆ the lack of enough measurements, any of these are violated, ˆ we take αLS and βLS to be zero. This means that, in this case, ˆ we can not estimate the correlated part of the channel. Once the underlying parameters of our model are estimated, channel at position q ∈ K can be estimated as follows. We have the following for the probability distribution of ΥdB (q), conditioned on all the gathered measurements and the underlying parameters: f (ΥdB (q)|YQ , θ, α, β, σ 2 ) ∼ 2 ˜ N ΥdB,Q (q), σdB,Q (q) with ˜ ΥdB,Q (q) E ΥdB (q) YQ , θ, α, β, σ 2

968

IEEE TRANSACTIONS ON WIRELESS COMMUNICATIONS, VOL. 11, NO. 3, MARCH 2012
−20 −30 −40 −50 indoor channel prediction outdoor channel prediction

−1 = hT (q)θ + φT (q)Rtot,Q YQ − HQ θ and Q 2 σdB,Q (q)

E

˜ ΥdB (q) − ΥdB,Q (q)

2

θ, α, β, σ

2

−1 = α + σ 2 − φT (q)Rtot,Q φQ (q), Q

(7) q − qb
ANMSE (dBm)

−60

where h(q) = 1

, · · · , e− denotes the cross and φQ (q) = α e− covariance between Q and q. Therefore, the Minimum Mean Square Error (MMSE) estimation of ΥdB (q), assuming perfect ˜ estimation of the underlying parameters, is given by ΥdB,Q (q). We then have the following by considering the true estimated ˆ ˆ ˆ ˆ −1 ˆ parameters: ΥdB,Q (q) = hT (q)θ + φT (q)Rtot,Q YQ − HQ θ , Q
T ˆ ˆ ˆ ˆ ˆ and where φQ (q) = α e− q−q1 /β , · · · , α e− q−qk /β ˆ ˆ tot,Q = αRnorm,Q (β) + σ 2 Ik×k . R ˆ ˆ The prediction quality at position q improves, the more correlated the available channel measurements become with the value of the channel at position q. In order to mathematically assess this, the next lemma characterizes the average number of the available measurements at the β neighborhood of the point to be predicted, for the case of randomlydistributed available channel measurements in 1D. The βneighborhood of a point q, in the workspace K, is defined as {z ∈ K|d(z, q) < β}, where d(z, q) denotes the Euclidian distance between points z and q. Lemma 1: Consider the case that k channel measurements, at positions {q1 , q2 , · · · , qk } are available, for predicting the channel at point q. Let Nβ (Q, q) represent the number of points in Q = {qi }k , which are located in the β neighbori=1 hood of q, where q and {qi }k are i.i.d. random positions, i=1 uniformly distributed over the workspace K = [0, L]. We β β2 then have, Nβ (Q, q) = k 2 L − L2 , where Nβ (Q, q) = T

q1 −q β

− D{q}

T

, D{q} = 10 log10 qk −q β

−19 −20 −21 −22 −23 non adaptive (R1−R2) adaptive (R1−R2)

−22 −24 −26 −28 −30 5 10 15 20 25 30 non−adaptive (R2−R3) adaptive (R2−R3)

% of measurements

Fig. 3: Impact of different environments on channel prediction performance, using real channel measurements. (top) indoor and outdoor, (middle) main room (R1) and hallway (R2) of Fig. 10 and (bottom) hallways R2 and R3 of Fig. 10.

EQ,q Nβ (Q, q) and EQ,q {.} represents the expected value w.r.t. Q and q. Proof: The proof is straightforward. Special case - probabilistic path loss: If the knowledge of beta is not available or is not used in the prediction (thus beta is assumed zero), then Eq. 7 results in the same probability distribution for all the points that are equallyspaced from the transmitter. An example of this case can be seen in Fig. 1 (top-middle), where we have the same predicted probability of connectivity (probability that the ΥdB (q) is above a given threshold) for all the points at a given radius from the transmitter. Our more general case of Eq. 7 is then shown in Fig. 1 (top-right), where a probability distribution (and a resulting probability of connectivity) is assigned to each point in the workspace. Both these cases result in a more comprehensive channel prediction than the commonlyused disk model of Fig. 1 (top-left). Next we show the reconstruction of two real channels, using our proposed method. The performance metric is the Average Normalized Mean Square Error (ANMSE) of the estimated channel, where the following Normalized Mean Square Error, NMSE =
K

a street in downtown San Francisco [19] as well as for an indoor channel measurement along a route in the basement of the ECE building at UNM. The indoor experiment uses an 802.11g WLAN card while the outdoor measurement is based on measuring receptions from an AT&T cell tower [19]. For both cases, all the underlying parameters are estimated using the LS approach of this section. Consider the outdoor case, for instance. We have the measurements of the received signal power, every 2mm along a street of length 16m in San Francisco, mounting to 8000 samples. Fig. 3 (top) then shows the prediction performance where only a percentage of the total samples were available to a node. The available measurements are randomly chosen over the street. 5% measurements, for instance, means that a robot has collected 400 samples, randomly over that street, based on which it will predict the channel over the whole street. The prediction error variance is −29dB for the case of 5% measurements. It can be seen that both channels can be reconstructed with a good quality. The outdoor channel, however, can be reconstructed with a considerably better quality. This is expected as the indoor channel suffers from a more severe multipath fading, which makes it less spatially predictable. A. Space-varying Underlying Parameters and Adaptive Channel Prediction So far, we considered channel prediction over a small enough space such that the underlying channel parameters can be considered constant over the workspace. However, if the available channel samples belong to a large enough space (such as the entire floor), the underlying parameters can be space-varying. In this part, we show how the previous framework can be extended to an adaptive approach, in order to address the case where the operation, and the corresponding available channel measurements, are over a large space. Basically, a robot can use its localization and mapping information (which it will have for navigation and collision avoidance) to detect when something changes in the structure of its environment. For instance, it can detect when

ˆ (ΥdB (q)−ΥdB,Q (q))2 dA Υ2 (q)dA dB K

, is averaged

over several different randomly-selected sampling positions, for a given percentage of collected samples. Fig. 3 (top) shows the reconstruction performance for an outdoor channel across

MALMIRCHEGINI and MOSTOFI: ON THE SPATIAL PREDICTABILITY OF COMMUNICATION CHANNELS

969

NMSE of nPL estimation

it moves out of a room to a hallway or when it reaches an intersection. Thus, we assume that the underlying parameters can possibly change when some environmental factors change. From analyzing several real measurements, this is a reasonable assumption. While there could possibly be cases that are not captured by this assumption, i.e. having a drastic change in an underlying parameter without any environmental change, such cases are rare and the robot can not know about it to adapt its strategy anyways. Let a region denote a place of operation where there is no environmental changes and the underlying parameters can be considered constant (such as a room or a hallway with no intersection that leads to the transmitter). In order to allow the node to give less weight to the available measurements that are collected in different regions and/or are far from the position where the channel needs to be estimated, we introduce a forgetting factor and a distance-dependent weight. This allows the node to adapt the impact of a sample measurement on its prediction framework. The forgetting factor is used to let the node give less impact to a measurement if it belongs to a different region, as compared to the place where the robot needs to predict the channel. On the other hand, the distance-dependent weight allows the robot to give less weight to the farther measurements. Consider the case where the workspace consists of p different regions, i.e. K = p Ri . i=1 Let τi represent the region, where the ith measurement belongs to, i.e. qi ∈ Rτi . Define the forgetting matrix F , with the following characteristics: 1) F is symmetric, 2) F is stochastic and 3) [F ]i,j = fi,j is proportional to the similarity between regions i and j. The third property implies that, maxj fj,i = fi,i and fi,j ≥ fi,k iff regions i and j have more environmental features in common, as compared to regions i and k. Furthermore, let G denote the functional space of all non-increasing functions on R+ . For q ∈ Rm ⊂ K, we define the corresponding weight matrix as: [ΨQ (q)]i,i = fτi ,m × gτi ,m ( q − qi ) and [ΨQ (q)]i,j = 0 for i = j, where gτi ,m ∈ G. One candidate for g is an exponential function: gτi ,m ( q−qi ) = e bτi ,m . fτi ,m and bτi ,m are design parameters, which the robot can choose. They impact how conservative the robot will be in taking the measurements of different 1 2 ˆ regions into account. Let θWLS (q) = minθ ΨQ (q) YQ − HQ θ denote the weighted LS estimation of the path loss parameters, for prediction at position q ∈ Rm . We then have, −1 T T ˆ HQ ΨQ (q)YQ . The channel and θWLS (q) = HQ ΨQ (q)HQ other underlying parameters can be similarly estimated. Fig. 4 shows the performance of our adaptive approach when a robot moves along a street. The channel measurement is in reception from an AT&T cell tower, in a street in San Francisco [19], which experiences very different path loss exponents due to the presence of an intersection that leads to the transmitter. The robot samples the channel as it moves along the street and estimates the path loss slope, without any a priori information in this environment. The figure compares the performance of the non-adaptive case with that of the adaptive one and shows that we can benefit considerably from the adaptation. Next, Fig. 3 (middle) shows the prediction quality when a number of robots operate in our basement,
2 − q−qi Received power (dBm)

−50 −60 −70 −80 −90 10 1 10 0 10
2

log of distance (dB)

10

−5

non−adaptive adaptive
2.38 2.39 2.4 2.41 2.42 2.43 2.44

2.37

log of distance (dB)

Fig. 4: Performance of our adaptive approach, in estimating the path loss slope, when a robot moves along a street in San Francisco and samples the channel along its trajectory [19]–(top) channel received power across the street along with its best slope fit and (bottom) prediction error variance of the robot, as it moves along the street and measures the channel. over a large area and cooperate for channel prediction. The regions of operation are R1 and R2, as indicated in Fig. 10. Note that the performance is simulated, in this case, using real channel measurements in this environment. It can be seen that the adaptive approach can improve the performance as compared to the non-adaptive case. In the non-adaptive case, all the gathered and communicated measurements are utilized by each robot for channel prediction, without taking into account that these measurements may belong to different regions. It can be seen that we can benefit a couple of dBs, by using the adaptive approach. In other tests in different environments, we also observed that the adaptation may make a negligible difference if different regions are not that much different, in terms of their underlying parameters, as expected. Fig. 3 (bottom) shows an example of such a case for operation over a different area in our basement. It can be seen that the performance curves are very close. In this paper, it is our goal to fundamentally understand the impact of different environments (in terms of their underlying parameters) on the proposed channel prediction framework. Consider the four marked regions of Fig. 10 for instance. We want to understand how the channel prediction quality changes (and justify the observed behaviors) when we move from one region to another. Therefore, in the rest of the paper, we consider the non-adaptive channel prediction framework, to predict the channel over a region where the underlying parameters can be considered constant. We then characterize the impact of different environments (in terms of the underlying channel parameters) on the prediction framework. III. I MPACT OF C HANNEL PARAMETERS ON THE P REDICTION E RROR VARIANCE In this section, we characterize the impact of the underlying channel parameters on the spatial predictability of a wireless channel. We assume that the underlying parameters are estimated perfectly in this section to avoid error propagation from parameter estimation to channel prediction. In the subsequent sections, we then extend our analysis to take the impact of the estimation error of key underlying parameters into account.

970

IEEE TRANSACTIONS ON WIRELESS COMMUNICATIONS, VOL. 11, NO. 3, MARCH 2012

Let ΥdB (q) = 10 log10 Υ(q) represent the received signal strength at position q ∈ K in dB. Based on the gathered measurements at Q ⊂ K, the goal is to estimate the channel ˜ at q ∈ K\Q, using the channel predictor, ΥdB,Q (q) of Eq. 7, 2 with the corresponding error covariance of σdB,Q (q). We next characterize the impact of different channel parameters on this prediction. We first introduce the following lemmas. Lemma 2: Let Ψ(t) be an invertible matrix for t ∈ R. We −n n have dΨ = −Ψ−n dΨ Ψ−n , where n is a positive integer. dt dt Proof: Taking the derivative from both sides of equation Ψn (t)Ψ−n (t) = Ik×k , with respect to t, proves the lemma. Lemma 3: Let J be an n-by-m matrix with the rank of m and Ψ be an n-by-n full rank matrix. If matrix Ψ is positive definite (Ψ 0), then J T ΨJ is positive definite. Proof: See [20] for a proof. 2 Theorem 1: The estimation error variance, σdB,Q , is an 2 2 increasing function of α and σ for α, σ ∈ [0, ∞) and an invertible Rnorm,Q . Proof: We first show that the estimation error variance is an increasing function of σ 2 . Let φnorm,Q (q) = 1 α φQ (q) denote the normalized cross covariance be2 tween Q and q. We have σdB,Q (q) = α + σ 2 − −1 σ2 T αφnorm,Q (q) Rnorm,Q + α Ik×k φnorm,Q (q). For α = 0, we have α = 0, taking the derivative with respect to σ 2 (using Lemma d 2 2) and then applying Lemma 3 result in: dσ2 σdB,Q (q) = 2 −2 σ φnorm,Q (q) > 0, ∀σ 2 ∈ 1 + φT norm,Q (q) Rnorm,Q + α Ik×k [0, ∞) and for an invertible Rnorm,Q , which completes the 2 proof. We next prove that σdB,Q (q) is an increasing function of 2 α. First assume that σ = 0. Taking the derivative with respect d 2 to α results in: dα σdB,Q (q) = 1 − φT norm,Q (q) Rnorm,Q + d 2 dσ2 σdB,Q (q) α=0

function of α in this case too. We next characterize the impact of β on the prediction quality, using properties of the Euclidean Distance Matrix (EDM) [22]. Given the position set Q = {q1 , q2 , · · · , qk } ⊂ K, the EDM Π = πi,j ∈ Rk×k is defined entry-wise as Π i,j = πi,j = qi − qj 2 for i, j = 1, 2, · · · , k. We have the following properties for the EDM: √ √ 1) πi,j ≥ 0 for i = j and πi,j = 0 for i = j. √ √ 2) πi,j = πj,i . √ √ √ 3) πi,l + πl,j ≥ πi,j for i = j = l. Theorem 2: Matrix Π = πi,j ∈ Rk×k is EDM if and only if −VkT ΠVk 0, ΠT = Π and πi,i = 0 for 1 ≤ i ≤ k, where Vk is the full-rank skinny Schoenberg auxiliary matrix: −1T 1 k−1 √ Vk ∈ Rk×k−1 . 2 I(k−1)×(k−1) Proof: Readers are referred to [22] for the details of the proof. Theorem 3: Let T = ti,j ∈ Rk×k represent the entry2 wise square root of Π = πi,j ∈ Rk×k where ti,j = πi,j . If Π is EDM, then T is EDM. This case is of interest because it corresponds to the absolute distance matrix. Proof: Readers are referred to [22]–[24] for the details of the proof. Lemma 4: The Hadamard product (Schur product) of two positive-definite matrices is positive-definite and the Hadamard product of two positive-semidefinite matrices is positive-semidefinite. Proof: Readers are referred to Theorem 7.5.3 of [20] for more details. Theorem 4: The estimation error variance is a decreasing function of β ∈ (0, ∞) for σ 2 = 0 and a non-increasing function of β ∈ (0, ∞) for σ 2 = 0 and an invertible Rnorm,Q . Proof: Case of α = 0 is not of interest in this theorem since we are interested in the impact of shadowing. Therefore, in this proof we assume that α = 0. Let T q1 − q , q2 − q , · · · , qk − q represent δQ (q) = the distance vector between the set Q and position q ∈ Q diag δQ (q) . Let TQ i,j = qi − qj , and ΔQ (q) ∀qi , qj ∈ Q, denote the absolute distance matrix corresponding to the set Q. First assume that σ 2 = 0. We have Eq. 8, shown at the top of the next page, where (•) denotes the Hadamard product. Moreover, it can be conT firmed that Rnorm,Q ΔQ (q) = Rnorm,Q • 1k δQ (q) . There−1 d 2 1 T fore, we have: dβ σdB,Q (q) = − β 2 φT (q)Rtot,Q α 1k δQ (q) + Q 1

= 1 > 0, ∀σ 2 ∈ [0, ∞). For

2 2 −1 −2 σ2 φnorm,Q (q). Define + σ Rnorm,Q + σ Ik×k α Ik×k α α d 2 σdB,Q (q). f is of class C ∞ on R+ with the f (α) dα

we By using these three properties, we have f (α) > 0, which means that the estimation error variance is an increasing function of α ∈ [0, ∞). Furthermore, if σ 2 = 0 and Rnorm,Q is invertible, d 2 then dα σdB,Q (q) σ2 =0 = f (∞), which is positive as shown for property 2. Therefore, estimation variance is an increasing

−2 σ2 σ2 d φnorm,Q (q) = −φT norm,Q (q) α dα Rnorm,Q α Ik×k −2 σ2 φnorm,Q (q) = −2σ 4 φT norm,Q (q) αRnorm,Q α Ik×k −3 2 φnorm,Q (q). Since αRnorm,Q + σ 2 Ik×k σ Ik×k 0, d can then easily show that dα f (α) < 0 using Lemma 3.

following properties: 1) f (0) = 1, 2) f (∞) = 1 − −1 d φT norm,Q (q)Rnorm,Q φnorm,Q (q) > 0 and 3) dα f (α) < 0. First property can be easily confirmed. We next prove the second property. Let Rnorm,Q {q} represent the correlation matrix corresponding to Q {q}. We have Rnorm,Q {q} = Rnorm,Q φnorm,Q (q) , where Rnorm,Q {q} is assumed 1 φT norm,Q (q) invertible. Thus, under the assumption that Rnorm,Q is invertible, the second property can be easily confirmed, using the Schur complement of Rnorm,Q block [21]. d Next we prove the third property. We have dα f (α) = 2 2 −1 d σ + σ Rnorm,Q + −φT norm,Q (q) dα Rnorm,Q + α Ik×k α + +

−1 δQ (q)1T − TQ • Rnorm,Q + 2σ 2 ΔQ (q) Rtot,Q φQ (q). From k Lemma 4, we know that the Hadamard product of two positive-semidefinite matrices is positive-semidefinite. Thered 2 fore, to prove that dβ σdB,Q (q) σ2 =0 < 0, it suffices to show T that 1k δQ (q) + δQ (q)1T − TQ is positive-semidefinite we k T 0 δQ (q) 0 . Let T{q} Q = ∈ know that ΔQ (q) δQ (q) TQ R(k+1)×(k+1) represent the distance matrix corresponding to {q} Q. Let ei denote a unit vector in Rk+1 , where all the entries are zero except for the ith one. Therefore, the Schoenberg auxiliary matrix can be represented as Vk+1 = 1 √ e2 − e1 , · · · , ek+1 − e1 . We have: 2

MALMIRCHEGINI and MOSTOFI: ON THE SPATIAL PREDICTABILITY OF COMMUNICATION CHANNELS

971

d 2 1 −1 −1 −1 −1 σdB,Q (q) = − 2 φT (q) ΔQ (q)Rtot,Q − αRtot,Q TQ • Rnorm,Q Rtot,Q + Rtot,Q ΔQ (q) φQ (q) dβ β Q 1 −1 −1 = − 2 φT (q)Rtot,Q Rtot,Q ΔQ (q) − αTQ • Rnorm,Q + ΔQ (q)Rtot,Q Rtot,Q φQ (q) β Q 1 −1 −1 = − 2 φT (q)Rtot,Q α Rnorm,Q ΔQ (q) + ΔQ (q)Rnorm,Q − TQ • Rnorm,Q + 2σ 2 ΔQ (q) Rtot,Q φQ (q), β Q

(8)

d T −Cθ,ML dα HQ αRnorm,Q + σ 2 Ik×k

−1

HQ Cθ,ML

=

−

T Vk+1 T{q}

Q Vk+1

1 T ei+1 − e1 T{q} Q ej+1 − e1 2 1 = − eT T{q} Q ej+1 − eT T{q} Q ej+1 − eT T{q} 1 i+1 2 i+1 1 = qj − q + qi − q − qi − qj 2 1 T 1k δQ (q) + δQ (q)1T − TQ . = k 2 i,j =−

i,j

Q e1

Then, matrix T{q} Q is EDM using Theorem 3. Therefore, T applying Theorem 2 for EDM T{q} Q results in: 1k δQ (q) + T T δQ (q)1k −TQ = −2Vk+1 T{q} Q Vk+1 0, which completes the proof. Next consider the case where σ 2 = 0. A similar d 2 ≤ 0, under the derivation will result in dβ σdB,Q (q) σ2 =0 assumption that Rnorm,Q is invertible. Therefore, the estimation error variance is a non-increasing function of β in this case. Note that path loss parameters, KdB and nPL , do not affect the estimation error variance in this case. In Section VI, we show the impact of different environments (with different underlying parameters) on channel predictability, using several measurements in our building. We next characterize the impact of the underlying parameters on the estimation of path loss parameters. IV. I MPACT OF C HANNEL PARAMETERS ON PATH L OSS E STIMATION In this section, we explore the effect of the underlying channel parameters on the estimation of path loss parameters. To provide a benchmark, we first consider the ML estimator of Eq. 3, where we assume that α, β and σ 2 are perfectly known. We then consider the Least Square estimator of Eq. 5 for a more realistic case, where α, β and σ 2 are not known at the time of estimating path loss parameters. Let ˆ ˆ θML = [KdB,ML nPL,ML]T denote the ML estimation of path ˆ loss parameters as denoted by Eq. 3. We have the following ˆ ˆ T = error covariance matrix: Cθ,ML = E θ−θML θ−θML
2 = [Cθ,ML ]1,1 and σnPL,ML = ˆ ˆ [Cθ,ML ]2,2 denote the ML estimation error variance of KdB,ML and nPL,ML respectively. We have the following Theorem. ˆ 2 2 Theorem 5: Both σK ˆ dB,ML and σnPL,ML are increasing funcˆ 2 2 tions of α and σ for α, σ ∈ [0, ∞) and an invertible Rnorm,Q . dB,ML θ,ML d = dα eT Cθ,ML e2 = eT dα e2 > 0. To show that 2 2 the estimation error of path loss parameters is an increasing dCθ,ML function of σ 2 , it suffices to show that dσ2 0. We have, −1 dCθ,ML T d 2 HQ Cθ,ML = dσ2 = −Cθ,ML HQ dσ2 αRnorm,Q + σ Ik×k −2 T HQ Cθ,ML 0, for Cθ,ML HQ αRnorm,Q + σ 2 Ik×k α, σ 2 ∈ [0, ∞) and an invertible Rnorm,Q . In general, the estimation error variance of path loss parameters does not have monotonic behavior as a function of β. To get a better understanding of the impact of correlation distance on the estimation of path loss parameters, we consider two extreme cases of β = 0 and β = ∞. More specifically, we characterize the optimum positions of the measurement points at both extremes and find the minimum achievable estimation error variance. A. Case of β = 0: In this case, Rtot,Q β = 0 = (α + σ 2 )Ik×k and the error covariance matrix of path loss parameters can be characterized as:

+ σ Ik×k Rnorm,Q αRnorm,Q + By using Lemma 3 and the assumption dCθ,ML that Rnorm,Q 0, we can easily see that dα 0. Let T T e1 = [1 0] and e2 = [0 1] denote unit vectors in R2 . We dCθ,ML d 2 d = dα eT Cθ,ML e1 = eT dα e1 > 0 and have: dα σK 1 1 ˆ
T Cθ,ML HQ αRnorm,Q −1 HQ Cθ,ML. σ 2 Ik×k d 2 ˆ dα σnPL,ML dB,ML 2

−1

dC

β→0

T lim Cθ,ML = (α + σ 2 )(HQ HQ )−1

= (α + σ 2 ) =

T α + σ2 DQ DQ T T DQ Ak DQ 1k DQ

k −1T DQ k

−1T DQ k T DQ DQ

−1

1 T DQ k , k

(9)

T −1 HQ Rtot,Q HQ

−1

2 , where σK ˆ

where Ak = kIk×k − 1k 1T . As can be seen, the estimation k error variances of both KdB and nPL are functions of sampling positions (Q). Lemma 5: Matrix Ak = kIk×k − 1k 1T has 0 and k as k eigenvalues with the multiplicity of 1 and k − 1 respectively. Let v1 ∈ span{1k } and v2 ∈ 1⊥ , where 1⊥ = v v T 1k = 0 . k k We have Ak v1 = 0 and Ak v2 = kv2 . Proof: The proof is straightforward and is omitted. 1⊥ 1k k Theorem 6: Let DQ and DQ denote the projection of ⊥ DQ to span{1k } and 1k subspaces respectively. The optimum 2 2 positioning, which minimizes both σK ˆ dB,ML and σnPL,ML for the ˆ case of β = 0, is
1k QPL,β=0 = arg max ||DQ ||2 , s.t. Q ⊂ K and DQ = 0. 2 opt Q

Proof: σ 2 Ik×k
−1

We HQ
−1

have

Cθ,ML

= =

T HQ αRnorm,Q + 1 α RQ . Taking dCθ,ML = in: dα

(10) Proof: We have the following optimum positioning in order to minimize the estimation error variance of

, where Rnorm,Q

the derivative with respect to α results

972

IEEE TRANSACTIONS ON WIRELESS COMMUNICATIONS, VOL. 11, NO. 3, MARCH 2012

KdB , using Rayleigh-Ritz theorem [20]: Qopt σ2
2 arg min s.t. Q⊂K σK ˆ dB,ML,β=0 = arg max s.t.

1k Q Q ⊂ K and DQ = 0 . This optimization problem can have multiple solutions, depending on the structure of the space, all of which achieve the minimum error vari2 ance of α+σ . Similarly, we have the following to minik mize the estimation error variance of nPL : Qopt = σ2 2 T arg min s.t. Q⊂K σnPL,ML,β=0 = arg max s.t. Q⊂K DQ Ak DQ = ˆ 2 arg max s.t. Q⊂K and D1k =0 ||DQ ||2 . Therefore, Eq. 10 repreQ sents the optimum positioning which satisfies both objectives. nPL,ML ,β=0 ˆ

ˆ KdB,ML ,β=0 T DQ Ak DQ Q⊂K DT DQ Q

= =

Next, we provide an intuitive interpretation. Similar to Eq. 2, the measurement vector can be represented by YQ = 1 (KdB × k)u1 + (−nPL DQ 2 )u2 + ΞQ , where u1 = √k k DQ and u2 = DQ 2 are normalized vectors. Then, the problem becomes similar to the decoding problem in CDMA (Code Division Multiple Access) systems. Thus, we have DQ ∈ 1⊥ . k Moreover, maximizing k and DQ 2 , which can be interpreted as maximizing the SNR of each term, results in a better estimation of KdB and nPL respectively. B. Case of β = ∞: Next we characterize the impact of correlation on the estimation quality of path loss parameters, when β goes to ∞. To simplify the derivations, we define α two variables: ρ = σ2 for σ 2 = 0, which denotes the ratio of the power of shadowing to multipath power (in dB) and χ = α + σ 2 , which represents the sum of the two powers. The following can be easily confirmed for σ 2 = 0: 1) χ −1 limβ→∞ Rtot,Q = ρ1k 1T + Ik×k 1+ρ , 2) limβ→∞ Rtot,Q = k ρ T 1+ρ Ik×k − 1+ρk 1k 1k χ (using Matrix Inversion Lemma), 3) 1+ρ 1 −1 −1 limβ→∞ 1T Rtot,Q 1k = k 1+ρk χ , 4) limβ→∞ 1T Rtot,Q DQ = k k
1+ρ 1 T = DQ ρAk + Ik×k DQ 1+ρk χ . Using the above equations, we have Eq. 11, shown at the top of next page. Remark 1: It can be seen from Eq. 11 that Theorem 6 also characterizes the optimum positioning for this case. Moreover, opt if QPL denotes the solution of Eq. 10, then we have, 2 σK ˆ dB,ML β=0 2 σnPL,ML ˆ 1+ρ 1 1+ρk χ , 5) 2 1+ρ ρ T 1+ρk 1k DQ χ

1 T DQ k

T −1 limβ→∞ DQ Rtot,Q DQ =

T DQ DQ −

distribution. Thus, for k ≥ 2, the slope of path loss, −nPL , can be perfectly estimated. However, the uncertainty of results in a bias in the estimation of KdB , as can be seen from Eq. 13. It can also be seen that the estimation error covariance is not a function of the sampling positions anymore. We next characterize the LS estimation of path loss parameˆ ters. Let θLS denote the LS estimation of path loss parameters as denoted by Eq. 5. We have the following error covariance −1 T −1 T T HQ Rtot,Q HQ HQ HQ . The matrix: Cθ,LS = HQ HQ following Theorem characterizes some properties of this estimator. ˆ Theorem 7: Let θLS and Cθ,LS represent the Least Square estimator of path loss parameters and the corresponding es2 timation error covariance matrix respectively. Let σK ˆ dB,LS and 2 ˆ dB,LS and nPL,LS σnPL,LS denote the LS error variances of K ˆ ˆ respectively. We have the following properties: 1) Cθ,LS Cθ,ML . 2 2 2 2) σK ˆ dB,LS and σnPL,LS are increasing functions of σ for α, ˆ 2 2 σ 2 ∈ [0, ∞). Moreover, σK ˆ dB,LS and σnPL,LS are increasing ˆ 2 functions of α for α, σ ∈ [0, ∞) and an invertible Rnorm,Q . 3) Both ML and LS estimators provide the same estimation error covariance matrices if β = 0 or ∞. Proof: The first property says that the ML estimator outperforms the LS one, as expected. We skip the mathematical proof due to space limitations. The second property can be easily confirmed by taking the derivatives with respect to σ 2 and α. We next prove the third property. For β = 0, we have Rtot,Q = (α + σ 2 )Ik×k , resulting in limβ→0 Cθ,LS = T DQ DQ 1 T DQ α+σ2 k . For β = limβ→0 Cθ,ML = DT A DQ T 1 k DQ k Q k T 2 ∞, we have Rtot,Q = α1k 1k + σ Ik×k , Cθ,LS =
T T T T H Q 1 k 1 T H Q HQ H Q α HQ HQ + σ 2 H Q HQ k and Eq. 14, shown at the next page. Therefore, we have T σ2 1 0 DQ DQ 1 T DQ k lim Cθ,LS = α + T T 0 0 k β→∞ DQ Ak DQ 1k DQ ⎤ ⎡ T T D DQ 1 D α + DT Q DQ σ 2 DTkA Q Q σ 2 Q Ak Q kD ⎦. =⎣ T 1k DQ k σ2 σ2 DT A DQ DT A DQ
Q k Q k

−1

−1

−1

β=0

α + σ2 2 σ2 , σK , = =α+ ˆ dB,ML k k β=∞ α + σ2 σ2 2 = and σnPL,ML = ˆ 2 DQopt DQopt β=∞
PL PL

2

. (12)

As can be seen, the fully correlated case provides a smaller estimation error variance for nPL and larger for KdB . In [5], we showed that the slope of path loss, nPL , has the most impact on the overall channel estimation error variance. Thus, case of β = ∞ would be more desirable than β = 0. Remark 2: Consider the case where multipath effect is negligible, i.e., σ 2 = 0. We have β→∞ lim Cθ,ML,σ2 =0 =

α 0

0 . 0

(13)

For this case, the measurement vector becomes YQ = HQ θ + 1k , where ∼ N (0, α) with N denoting a Gaussian

By comparing this equation to Eq. 11, the third property can be verified. Remark 3: Theorem 7 (3) shows that the optimum positioning of Eq. 10 minimizes the estimation error variance of the LS case too. We next verify the derived theorems, using a simulated channel. Fig. 5 shows a simulated channel, generated with our probabilistic channel simulator [25], with the following parameters: frequency of operation of 1GHz, θ = [−22 3.0]T , √ √ α = 8dB and β = 1m. As for multipath fading, this channel experiences a correlated Rician fading, with Jakes power spectrum [8], which results in the multipath fading getting uncorrelated after 0.12m. The pdf of a unit-average Rician distribution, with parameter Kric , is given by [6]: fΥMP (x) = (1 + Kric )e−Kric −(1+Kric )x I0 2 xKric (Kric + 1) , where I0 (.) is the modified zeroth-order Bessel function. Note that Kric = 0 results in an exponential distribution, which

MALMIRCHEGINI and MOSTOFI: ON THE SPATIAL PREDICTABILITY OF COMMUNICATION CHANNELS

973

β→∞

lim Cθ,ML = lim =

β→∞

T −1 HQ Rtot,Q HQ

−1

= lim

β→∞

−1 1T Rtot,Q 1k k −1 −1T Rtot,Q DQ k T DQ

−1 −1T Rtot,Q DQ k −1 T DQ Rtot,Q DQ

−1

1 1 + ρk χ× × T 1+ρ 1 + ρk DQ Ak DQ ⎡ ⎤ DT DQ 1T D α + DT Q k DQ σ 2 DTkAkQ Q σ 2 D QA Q ⎦. =⎣ 1T DQ k 2 k σ σ2 T T D A DQ D A DQ
Q k Q k

ρAk + Ik×k DQ 1 T DQ k

1 T DQ k k

=

χ 1+ρ

⎡ ρ+ ⎣

T DQ DQ T DQ Ak DQ 1T DQ k T DQ Ak DQ

1T DQ k T DQ Ak DQ ⎦ k T DQ Ak DQ

⎤

(11)

T (HQ HQ

−1

T HQ 1 k 1 T HQ k 2 T DQ DQ T 1 k DQ

T HQ HQ

−1

=

1
T DQ Ak DQ

1 T DQ k k

k2 −k 1T DQ k
2

−k 1T DQ k 2 1 T DQ k

T DQ DQ T 1 k DQ

1 T DQ k k
3 T DQ DQ T 1 k DQ

=

1
T DQ Ak DQ 2

T k 2 DQ DQ − k 1 T DQ k 0

T −k DQ DQ 1T DQ + 1T DQ k k 0

1 T DQ 1 k = k 0

0 0

(14)

9
Received power (dBm)

+2/k=8.25 (Eq. 13)

30 20 10 −10 −20 −30 1
K
2

8 7 6 dB 0

5 4 3

random pos. with ML est. random pos. with LS est. opt. pos. of Theorem 6

(+ )/k=1.25 (Eq. 13)

2

0.5

0 y (m)

−0.5

−1 −1

−0.5

0 x (m)

0.5

1

2 1 −4 10 10
−2

10

0



10

2

10

4

following underlying parameters: θ = [−22 3.0]T , α = 8dB, √ β = 1m and σ = 2dB. The transmitter is located at qb = [0 0]T .

Fig. 5: A 2D simulated channel at 1GHz frequency with the √ √

Fig. 6: Impact of β on the estimation of KdB for both optimum positioning of Theorem 6 and random sampling.
0.35 0.3 0.25
PL

random pos. with ML est. random pos. with LS est. opt. pos. of Theorem 6
2 2 2 2  /||DQopt.|| =0.0091 (Eq. 13)
PL

experiences a considerable amount of channel variations, while Kric = ∞ results in no fading, i.e., we will have a channel with only path loss and shadowing. Multipath power (in dB), σ 2 , is related to Kric as follows: σ 2 = E ω 2 (q) = log10 (x)fΥMP (x)dx . √ For the simulated channel of Fig. 5, σ = 2 dB, which corresponds to Kric = 19. Fig. 6 and Fig. 7 show the impact of the correlation distance, β, on the estimation variance of KdB and nPL respectively. In this example, the workspace is a ring with an inner radius of 0.3m and an outer radius of 3.3m, superimposed on the simulated channel of Fig. 5, such that the centers of the rings are positioned at the transmitter. We consider the case where k = 8 samples are taken from the workspace. Furthermore, we compare the performance for the case of random uniformlydistributed samples with the case where samples are optimally positioned based on Theorem 6. For this workspace, enforcing 2 1k DQ = 0 results in max DQ 2 = 100k log2 ( 10 ), which can 10 3 100
∞ 0

2

0.2

n

(+ )/||DQopt.|| =0.0457 (Eq. 13)
PL

0.15 0.1 0.05 0 −4 10

log2 (x)fΥMP (x)dx−100 10

∞ 0

2

10

−2

10

0



10

2

10

4

Fig. 7: Impact of β on the estimation of nPL for both optimum positioning of Theorem 6 and random sampling.

be achieved if and only if half of the samples are distributed on the inner circle while the other half are on the outer one. Therefore, we assume that four samples are equally-spaced on the inner circle while the other four are equally-spaced on the

974

IEEE TRANSACTIONS ON WIRELESS COMMUNICATIONS, VOL. 11, NO. 3, MARCH 2012

−24

Sensitivity w.r.t.  Sensitivity w.r.t.  Sensitivity w.r.t. n

pl dB

ANMSE (dB)

−26

Sensitivity w.r.t. K

−28

−30