Visible light communications (VLC) based positioning systems will form an important part of the future generation wireless communication systems because they oﬀer higher accuracy for indoor positioning as compared to radio frequency based systems. In this paper, we have used non Bayesian statistical signal processing techniques for the hybrid time-of-arrival/received-signal-strength and hybrid time-diﬀerence- of-arrival/received-signal-strength based positioning. These hybrid measurements are combined with the following fusion algorithms: weighted least squares and the best linear unbiased estimator . These two fusion algorithms are compared in terms of the mean Euclidean error as a function of various parameters such as signal-to-noise ratio, transmitter arrangement and synchronization error. Even though the performance of the weighted least squares algorithm is better, the best linear unbiased estimator is still an attractive algorithm for systems that require a lower complexity. Index Terms Bestlinearunbiasedestimator,meanEuclideanerror,received-signal-strength,time-diﬀerence-of-arrival, time-of-arrival, visible light communications based positioning, weighted least squares I. Introduction Light-emitting-diode (LED) based visible light communications (VLC) has attracted a lot of attention in recent years [1]. This is primarily because such systems can provide illumination as well has high speed data rates to the users located indoors. Furthermore, such VLC systems can also help to position the receiver terminals indoors with suﬃcient accuracy as opposed to WiFi based The author may be reached at muhammad.bashir@fulbrightmail.org or salman.bashir@alumni.purdue.edu for any helpful comments or suggestions. 2 localization that suﬀers from interference and fading due to multipath [2]. The accurate positioning performance provided by VLC systems is very important for future applications in sixth generation wireless systems such as robot navigation, as well as location based services and advertisements among many others [1]. A. Literature review The focus of this paper is on utilizing statistical signal processing techniques in order to localize the receiver in an indoor visible light communications setting. Therefore, for the purpose of literature review, we have focused mainly on research that has approached visible light positioning from a signal processing perspective. The authors in [2] have proposed Epsilon, which is a visible light positioning system based on the received signal strength measurements. Epsilon uses binary frequency shift keying signaling and channel hopping in order to avoid collisions from the uncoordinated LED sources. The positioning accuracy is reported to be 0.4 meters, 0.7 meters and 0.8 meters for diﬀerent rooms/locations. The authors in [3] derive the Cramer-Rao Lower Bound (CRLB) for the time-of-arrival (TOA) based transmitter/receiver distance estimation for windowed sinusoidal signals. The CRLB is derived under the assumption that the transmitter/receivers are synchronized in time. The authors in [4] improve upon this work by deriving the CRLB’s for TOA, received-signal-strength (RSS), and the hybrid TOA-RSS systems, and they consider both the synchronous and asynchronous scenarios in their work. Moreover, the maximum likelihood estimators (MLE) are also proposed for the RSS and TOA based distance estimation. Furthermore, they also consider the eﬀect of the sampling time on the CRLB. The same authors devised the positioning algorithm in [5] by extending the work in [4]. They consider a direct positioning algorithm in which the position is estimated by maximizing a likelihood function, and a two-step positioning algorithm in which the maximum likelihood estimates of time-delay and RSS are combined using a nonlinear least squares algorithm. Additionally, they consider the positioning for synchronous as well as asynchronous systems and derive the CRLB curves for diﬀerent positioning scenarios. In their work, they conclude that the direct positioning yields more accurate results than the two-step system. However, the better accuracy of direct positioning comes with the extra cost of higher computational complexity. The focus of [6] is to consider the eﬀect of multipath on the TOA based positioning. They have modeled the delay proﬁle as a Gamma function, and analyzed the eﬀect of the shape parameter and 3 the scale parameter of the Gamma function on the Generalized Cramer-Rao Lower Bound (GCRLB). However, in most VLC scenarios, the non-line-of-sight (NLOS) component of the received signal is insigniﬁcant compared to the line-of-sight (LOS) signal [3], [4], and hence the application of [6] is limited to the situations where the NLOS portion of the signal is signiﬁcantly larger. The authors in [7] have presented an angel-of-arrival based positioning algorithm that employs least squares to achieve localization. They have also employed a Kalman ﬁlter assuming that the motion model in two dimensions is a Gaussian process. However, the average error of their localization algorithm is 0.2 meters on average, which may be considerable for applications that require more accuracy (for example robots trying to enter a narrow passageway etc.). In [8], the authors have proposed positioning based on one LED source and multiple receivers, in which the multiple receivers in a room cooperate during the process of localization. However, their algorithm is only feasible (achieves an accuracy of about 0.05 meters) when the receivers lie in certain circular patterns. Finally, in [9], the authors have proposed a machine learning based localization algorithm. In their paper, the ﬁrst task is to carry out the RSS based ﬁngerprinting followed by training of classiﬁers based on these ﬁngerprinting results via machine learning techniques. The output of the classiﬁers are fused using robust least squares algorithms, and a singular valued decomposition (SVD) based method is proposed to deal with the numerical stability problem in case the prediction matrix is singular. However, ﬁngerprinting and training of classiﬁers is a tedious and time consuming process, especially if the whole building is to be mapped for the purpose of ﬁngerprinting. Moreover, the SVD based method also incurs signiﬁcant computational complexity [9]. Regarding beam tracking in free-space optical communications, we refer the interested readers to [10]–[21] B. Contributions of this paper We primarily build on the work [5] in order to devise positioning algorithms based on visible light communication systems. The work presented in our paper is novel because we propose multiple optical sensors associated with a given receiver instead of positioning with just one sensor. In this sense, we are fusing more independent pieces of information (the observations for any given sensor are generated independently of the other sensors). This results in an increase in the Fisher information of the estimation problem, or, in other words, a reduction in the CRLB. This implies that the variance of the proposed estimators for a multiple sensor problem will be much smaller than a single sensor 4 case. Hence, the larger the number of the receiver sensors, the more accurately we can position the receiver. Secondly, based on the received-signal-strength, we weight the measurements obtained for each transmitter. We used two types of weighting schemes in this regard: i) An RSS based and ii) a binary weighting scheme for the sources that have the largest signal strength at the receiver. Moreover, we consider three types of LED transmitter arrangements in a room, and compare their performance by simulations. In order to “fuse” the measurements obtained for each transmitter/receiver pair, we use two algorithms. The measurements for the fusion algorithms are the hybrid maximum likelihood estimates of time-of-arrival (TOA), time-diﬀerence-of-arrival (TDOA) and received-signal-strength (RSS) corre- sponding to each sensor/transmitter pair. The ﬁrst of these fusion algorithms is the weightednonlinear least squares (WLS) that fuses the hybrid TOA-RSS and the hybrid TDOA-RSS measurements in order to localize the receiver. The second proposed algorithm “tracks” the position of the receiver in that it estimates the change in the position of the receiver. If the change in the position at the current instant of time is small from the position at the previous instant, the nonlinear equations in range and time can be linearized over that small change by the Taylor Series approximation. This results in a linear model for which we use a low complexity estimator, such as the best linear unbiased estimator (BLUE), which will fuse the maximum likelihood estimates of the time-delay at the receiver corresponding to the transmitters with the highest RSS. Therefore, the BLUE also fuses the hybrid TOA-RSS measurements in order to estimate the change in position. The process of fusion of TOA /TDOA and RSS estimates using the two algorithms described above is a similar idea to the two-step positioning system described in [5]. However, the fusion algorithms devised in this paper are diﬀerent than those discussed in [5], and yield a better mean Euclidean error performance as argued in the ﬁnal part of the paper. Moreover, we left out the direct positioning algorithm in our paper because of the higher complexity associated with it [5]. Furthermore, we analyze the eﬀect of the following factors on the positioning performance of both the proposed algorithms: 1) signal-to-noise ratio, 2) total number of sensors located on the receiver, 5 T1T2 T3 Fig. 1. An example of positioning via Triangulation. The transmitters are present at locations T1,T2andT3.The diameters of the circles represent the estimated time-delays corresponding to each transmitter. These estimated time-delays are based on the TOA or TDOA or RSS based measurements. For the ideal case when the time-delay estimates are exactly equal to the true value (shown by solid circles), the position where the three circles intersect (shown by black dot) represents the exact location of the receiver. However, because of noise, the estimates from transmitter T2are not accurate enough (represented by the dashed red circle), and the three circles do not intersect at one point causing an error in the position estimate (erroneous position is shown by the blue dot). 3) the total number of LED transmitters in a given location and their arrangement, 4) and the time synchronization error. C. Model assumptions For visible light based position estimation, we make the following assumptions for our model: 1) Each transmitter is assigned a unique frequency (or wavelength) band in which it can transmit. Similarly, the receiver can diﬀerentiate the signals from diﬀerent transmitters based on their wavelengths, and ﬁlter those into separate channels for processing. 2) A pilot channel is dedicated solely for the purpose of positioning. 3) A rectangular pulse is used in order to position the terminal. However, the discussion is general, and holds for any pulse shape. 4) It is assumed that the power from the non-line-of-sight components of the received light is negligible for the purpose of positioning [3], [4]. Hence, only the line-of-sight (LOS) component is considered for our positioning system. II. Problem formulation and measurement models for visible light positioning Let us assume that there are there are NTlight-emitting-diode (LED) transmitters located at points p(i) t,i= 1,2,...NT,and each p(i) t∈R3.Let us further assume that the receiver is made up of 6 Msensors, and those sensors are arranged symmetrically on the device (or a robot etc.) that needs to be tracked. If the location of sensor misp(m) r∈R3form= 1,...,M, then, roughly, the center of the device (the point that we want to track) is the sample average: pr,1 MM/summationdisplay m=1p(m) r, (1) where pris the location of the center of the receiver. If we assume that the height zof the object from the ground is ﬁxed, then we only want to track thexandycoordinates of its position on the ﬂoor. Then pm—and therefore pr—is just a function ofxandybecausezis a known ﬁxed constant. The received signal intensity (for the LOS scenario) at the mth sensor due to the ith transmitter is given by λ(im) r(t) =αimλ(i) s(t−τim) +λ(m) 0, (2) wheretinside some observation interval [0,T], αimis the attenuation factor of the optical channel, λ(i) s(t)is the transmitted signal intensity of the ith transmitter, and λ(im) r(t)is the received signal intensity at the mth receiver sensor due to the ith transmitter. The parameter λ(m) 0is a constant that represents the noise intensity at the mth receiver sensor due to thermal eﬀect and background radiation. The variable τimrepresents the TOA parameter [4], and it represents the delay introduced by theith transmitter mth sensor channel τim,T0+Rim c, (3) fori= 1,2,...,NT,m= 1,...,M. In (3),Rimrepresents the distance between the ith transmitter and themth receiver sensor, crepresents the speed of light in air in meters/sec, and T0is the time at which all the transmitters begin the transmission of the pulse. If the receiver is synchronized with the transmitters, then T0will be considered a known quantity. We now denote the channel attenuation factor and the time delay for the ith transmitter/ mth receiver pair by αimandτim, respectively. The quantity αimis modeled as [3] αim=/lscript+ 1 2πcos/lscript(φ) cos(θ)Ar R2 im, (4) where/lscriptis the Lambertian order, Aris the area of the photo-detector at the receiver and φand 7 LED φ θ optical sensor Fig. 2. Geometry of the transmitter (LED) and receiver (optical sensor) in the visible light positioning problem. θrepresent the irradiation and incidence angles, respectively. However, when the LED transmitter is looking downwards, and the receiver is pointing upwards—which is usually the case as shown in Figure 2— then φ=θ=⇒cos(φ) = cos(θ) =h/Rim, wherehrepresents the height of the LED transmitter with respect to the VLC receiver [4]. In that case [4], αim=/lscript+ 1 2πh/lscript+1ArR−/lscript−3 im. (5) A. Measurements for position estimation In this section, we assume that the positioning system uses a pilot channel dedicated for the purpose of positioning. Hence, a positioning pulse of width Tpis transmitted every Tseconds (or an integer multiple of Tseconds) in order to estimate the position of the receiver. Analyzing (2), we see that there are two unknowns in that equation: the delay τim,and the attenuation factor αim.Both these factors will depend on the distance of a given sensor from a particular transmitter. Hence, we can use the measurements based on τimandαimto our advantage as we estimate the position of a receiver. In order to estimate the position of the receiver, we ﬁrst form the parameter vector (with respect to theith transmitter) as θi,/bracketleftbigg αi1αi2···αiMτi1τi2···τiM/bracketrightbigg| .Furthermore, based on the estimate ˆθiof the vector θifor alli, we can fuse these estimates in order to estimate the coordinates of the receiver position in two dimensions. 8 A t→r(t) αimAA t→αimr(t−τim) τimTp T Fig. 3. The transmitted pilot pulse from the ith transmitter reaches the mth sensor with a delay of τimseconds, and channel attenuation factor αim. The time starts at time T0. 1) Measurement of time-delay with maximum likelihood estimation In this section, we use the maximum likelihood principle in order to generate the measurements ˆτimfori= 1,···,Nt, m= 1,···,M,in order to estimate the position of the receiver. In this regard, we employ photomultiplier tube receivers, which are commonly used for the detection of signals in optical communications. Furthermore, we can estimate τi(the average time it takes for the light signal from the ith transmit- tertoreachthecenterofthereceiver)basedon ˆτim.Since,themaximumlikelihoodestimatecommutes over a linear transformation, the maximum likelihood estimate ˆτiis simply ˆτi=1 M/summationtextM m=1ˆτim. The values ˆτiwill form the measurements/observations for the WLS or BLUE based fusion algo- rithm that we will discuss later in this paper. Let us assume that the mth sensor of the receiver receives the following positioning pulse by the ith transmitter after the conversion of the optical signal into electric charge [4]: xim(t) =αimRpr(t−τim) +wim(t), (6) wherer(t)is the transmitted optical signal and wim(t)is a zero mean additive white Gaussian noise (AWGN) process (with two-sided power spectral density N0/2) that takes into account the eﬀect of shot noise , thethermal noise and thebackground radiation . The factor Rpis theresponsivity of the photo detector [4]. The signal xim(t)is then ﬁltered by a low pass ﬁlter of frequency cut-oﬀ BHz, whereBis the eﬀective bandwidth of the pulse s(t).Let us denote Rpr(t)bys(t).The resulting ﬁltered signal is then discretized in time to obtain xim[n] =αims(n∆−τim) +wim[n]. (7) 9 In (7) ∆is the sampling interval, and it is set equal to 1/(2B)seconds, the Nyquist rate . Moreover, xim[n],x/prime im(n∆)andwim[n],w/prime im(n∆)wherex/prime im(t)is the ﬁltered signal, and w/prime im(t)is the ﬁltered noise bandlimited to BHz. It can be shown that wim[n]is a zero mean white Gaussian noise process with variance σ2=N0B[?]. Equation (7) can be further simpliﬁed as xim[n] =αims[n−Dim] +wim[n], n = 0,1,...,N−1, (8) whereDimis an integer, and represents the delay in discrete time, and N,T/∆is the observation interval in discrete time. We assume that the sampling period ∆is small enough such that τim≈ Dim∆. The signal xim[n]is now passed through a discrete-time matched ﬁlter of impulse response s∗[−n] to obtain yim[n] =xim[n]∗s∗[−n] =αims[n−Dim]∗s∗[−n] +wim[n]∗s∗[−n] (9) =αimrss[n−Dim] +rws[n], n = 0,1,...,N−1, (10) whererss[n]andrws[n]are the autocorrelation and cross-correlation functions, respectively, and ∗ representsthediscrete-timeconvolutionoperator.Furthermore, rws[n]isazeromeanGaussianprocess with variance E(r2 ws[n]) =σ2EswhereEs,∞/summationdisplay n=−∞s2[n].TheN×Ncovariance matrix of rws[n]is C,σ2 rss[0] rss[1]···rss[N−1] rss[−1] rss[0]···rss[N−1] ............ rss[−(N−1)]rss[−(N−2)]···rss[0] . Our objective is to compute the maximum likelihood estimates ˆDimusing (10), and based on these estimates, compute the maximum likelihood estimates ˆτimusing ˆτim= ∆ ˆDim.Furthermore, ˆτiis just a sample average of ˆτim’s. Finally, the maximum likelihood estimate of the range Rimis,ˆRim=cˆτim. The maximum likelihood estimate of Dimbased on the observations in (10) is deﬁned as ˆDim= arg max Dimp(yim|Dim) = arg max Dim1/radicalBig (2π)N|C|exp/parenleftbigg −1 2(yim−µim)|C−1(yim−µim)/parenrightbigg ,(11) 10 where yim,/bracketleftbigg yim[0]yim[1]···yim[N−1]/bracketrightbigg| ,and µim,/bracketleftbigg αimrss[0−Dim]αimrss[1−Dim]···αimrss[N−1−Dim]/bracketrightbigg| . Solving (11) further, the maximum likelihood estimate of Dimcan be written in a compact form [?]: ˆDim= arg max n(αimrss[n−Dim] +rws[n]), (12) and therefore, the estimate ˆDimis just the instant at which the output of the matched ﬁlter is maximized during the observation interval T. 2) Estimation of signal strength with maximum likelihood estimation For hybrid RSS-TOA based position estimation, we will need to estimate the received signal strength, or, the received signal attenuation factor αimalong with the time-delay estimates ˆDim. Based on the estimated values of the signal strength, we will use the best ηsources, where ηis an integer and η≥3. Thus we use three or more sources that suﬀer the smallest attenuation in order to localize the receiver through triangulation /multilateration . In order to estimate the attenuation factor αim,one option is to use a suﬃcient statistic in order to estimate αimthat does not depend on the knowledge of Dimorτim.One way of eliminating the dependence on Dimis to average the observation yimin (10). Thus, we will determine the maximum likelihood estimate of αimbased on a sample average of yimduring the observation interval (−T,T): ¯yim,1 2N−1N−1/summationdisplay n=−(N−1)yim[n] =1 2N−1N−1/summationdisplay n=−(N−1)αimrss[n−Dim] +1 2N−1N−1/summationdisplay n=−(N−1)rws[n](13) =αim 2N−1N−1/summationdisplay n=−(N−1)rss[n] +1 2N−1N−1/summationdisplay n=−(N−1)rws[n] =αim¯rss+ ¯rws (14) where ¯rssand¯rwsrepresent the sample averages for rss[n]andrws[n], respectively. The ﬁrst equality in (14) is true if Nis large enough since the sum of rss[n−Dim]is the same as the sum of rss[n]over (2N−1)terms when Nis large. It should also be noted that ¯rwsis a zero mean Gaussian random variable. It can be shown easily that the variance of ¯rwsis σ2 ws,E((¯rws)2) =1 (2N−1)2(2N−1)/summationdisplay i=1(2N−1)/summationdisplay j=1[Cws]ij (15) where [Cws]ijis the (i,j)element of the (2N−1)×(2N−1)covariance matrix Cwsof the vector 11 /bracketleftbigg rws[−(N−1)]rws[−(N−2)]···rws[0]rws[1]···rws[N−1]/bracketrightbigg| . Using (14), the maximum likelihood estimate of αimis ˆαim= arg max αim1/radicalBig 2πσ2 wsexp/parenleftBigg −1 2σ2 ws(¯yim−αim¯rss)2/parenrightBigg . (16) Taking the derivative of the natural logarithm of the function in (16) with respect to αim,and setting that equal to zero yields the maximum likelihood estimate: ˆαim=¯yim ¯rss.If, however, ˆαim<0,then we setˆαimequal to zero1. III. Fusion of measurements Now that we have generated the measurements ˆαimandˆτim,we are ready to fuse these measure- ments using the two algorithms discussed below in order to position the receiver. In the following analysis, we assume that the true position of the receiver is (xr,yr), and the estimate output by the fusion device is (ˆxr,ˆyr). A. Weighted least squares estimate of position We use the weighted least squares (WLS) approach for the nonlinear position estimation problem. In this regard, we ﬁrst deﬁne ˆαi,1 MM/summationdisplay m=1ˆαim. (17) Furthermore, we deﬁne the set Qwhich contains η(ηbeing a positive integer) largest ˆαi’s (corre- sponding to the ηbest sources). Thus, if/bracketleftbigg ˆα(1)ˆα(2)··· ˆα(NT)/bracketrightbigg| is the ordered array of attenuation factors ˆαisuch that ˆα(1)>ˆα(2)>···>ˆα(NT),thenQ,{ˆα(1),ˆα(2),..., ˆα(η)}where 0<η≤NT,and η∈Z+. In order to compute the WLS estimate, we minimize the function J(lr)with respect to vector lr: ˆlr= arg min lrJ(lr) = arg min lr(r−h(lr))|W(r−h(lr)), (18) where l(r)is the location vector of the receiver sensors deﬁned by lr,/bracketleftbigg p(1) rp(2) r...p(M) r/bracketrightbigg| ,and p(m) ris the location of the mth receiver sensor. The nonlinear function h(lr),which is an (Nt·M)×1 vector, is deﬁned as 1The value of αimequal to zero maximizes the likelihood function over the positive values of αim(theαim>0is the admissible range) since the likelihood function is quadratic in αim. 12 h(lr),/bracketleftBig ||p(1) t−p(1) r||2||p(2) t−p(1) r||2 . . .||p(NT) t−p(1) r||2||p(1) t−p(2) r||2 . . .||p(NT) t−p(NT) r||2/bracketrightBig| ,(19) where p(i) tis the position of the ith transmitter, and ||x||2represents the /lscript2norm of a given vector x. 1) Synchronous case Forthesynchronouscase,thequantity risdeﬁnedtobethe (Nt·M)×1vector/bracketleftbigg ˆR11ˆR21··· ˆRNt1ˆR12ˆR22··· ˆRNtM/bracketrightbigg| , where ˆRimis the estimate of the range between the ith transmitter and the mth sensor. The (Nt·M)×(Nt·M)diagonal matrix W(matrix of weights) is by deﬁnition W, ω110··· 0 0ω12··· 0 ............ 0 0···ωNtM , (20) whereωim,ˆαim/summationtextη i=1/summationtextM m=1ˆαimifˆαi∈Q,andωimis zero otherwise. In other words, we use the η“best” LED sources which suﬀer the smallest attenuation at the receiver, and then weight them according to their received signal strength. We can now rewrite (18) as ˆlr= arg min lrNT/summationdisplay i=1M/summationdisplay m=1/parenleftBigˆRim−||p(i) t−p(m) r||2/parenrightBig2ωim. (21) 2) Asynchronous case For the asynchronous case, we deﬁne the time-diﬀerence-of-arrival (TDOA) estimates as ˆtim,ˆτim−ˆτ(i−1)m, i= 2,3,...,NT,∀m. (22) The time delay estimate can be written as ˆτim=T0+τim+/epsilon1im,whereT0is the time at which the transmitters transmit the pulse, and /epsilon1imis the estimation error. Hence, ˆtim=T0+τim+/epsilon1im−(T0+τ(i−1)m+/epsilon1(i−1)m) =τim−τ(i−1)m+ (/epsilon1im−/epsilon1(i−1)m),(23) 13 1 2 iNTRiRni(xn,yn)(xn+1,yn+1) (ˆxn,ˆyn)(ˆxn+1,ˆyn+1) δxn+1δyn+1 · · ···· Fig. 4. Geometry of the visible light positioning problem for the best linear unbiased estimator. and the dependence on T0is eliminated. The weighted least squares estimate of the receiver position based on the TDOA is ˆlr= arg min lrNT/summationdisplay i=2M/summationdisplay m=1/parenleftBig/bracketleftBigˆRim−ˆR(i−1)m/bracketrightBig −/bracketleftBig ||p(i) t−p(m) r||2−||p(i−1) t−p(m) r||2/bracketrightBig/parenrightBig2ωim,(24) whereωim= 1if the corresponding ˆαiis inQ,and zero otherwise. Hence, we are using a binary weighting scheme for the asynchronous case. Finally, the estimate of the center of the receiver, ˆpr, for both the synchronous and asynchronous schemes, is given by ˆpr,1 MM/summationdisplay m=1ˆp(m) r. B. Best linear unbiased estimate (BLUE) of position (tracking) Using a linear model (linearity holds if the change in the position of the receiver is very small with respect to the sampling rate), we will obtain the position information as the best linear unbiased estimate (BLUE) of the true position based on the measurements ˆτi. The measurements ˆτican be written as ˆτi=T0+Ri/c+/epsilon1i, i= 1,2,...,NT, (25) where/epsilon1iis the measurement noise associated with ˆτi.It will be reasonable to assume that /epsilon1iis uncorrelated with /epsilon1jfori/negationslash=jandi,j= 1,2,...,NT.The variance of /epsilon1iis denoted by σ2 i.The quantityRirepresents the actual distance of the ith transmitter to the center of the receiver at the actual location (xn+1,yn+1)at timen+ 1. The geometry of this problem is depicted in Figure 4. 14 The nominal position (ˆxn,ˆyn)is the estimate of the true receiver location (xn,yn)at the discrete- timen. At time instant n+ 1, the distance Riof the receiver from the ith source is by deﬁnition Ri,/radicalBig (xn+1−xi)2+ (yn+1−yi)2.The Taylor expansion of Riabout the nominal distance Rni(the distance at instant n) gives Ri=Rni+xn−xi Rniδxn+1+yn−yi Rniδyn+1+ (higher order terms ). (26) Let us assume that the nominal position estimate (ˆxn,ˆyn)is accurate enough (close to (xn,yn)), and that the pulse transmission rate 1/Tis large enough compared to the movements of the receiver. Furthermore, let us assume that ˆRni≈Rni, that is, we have at our disposal a good estimate of the range Rni.Under this assumption, δxn+1andδyn+1will be small enough for linearization approximation to hold: Ri≈Rni+xn−xi Rniδxn+1+yn−yi Rniδyn+1. (27) If we assume that ˆRni≈Rniand(ˆxn,ˆyn)≈(xn,yn),we have that Ri≈ˆRni+ˆxn−xi ˆRniδxn+1+ˆyn−yi ˆRniδyn+1, (28) (29) and it follows that ˆRi=ˆRni+ˆxn−xi ˆRniδxn+1+ˆyn−yi ˆRniδyn+1+/epsilon1i. (30) (31) Therefore, in terms of time-delay ˆτi≈T0+ˆRni c+ˆxn−xi ˆRnicδxn+1+ˆyn−yi ˆRnicδyn+1+/epsilon1/prime i∀i, (32) where/epsilon1/prime i,/epsilon1i/c.For a synchronous system, the quantity T0is known. For this case, let us deﬁne ξi,ˆτi−T0−ˆRni/c,and the linear model is ξi=ˆxn−xi ˆRnicδxn+1+ˆyn−yi ˆRnicδyn+1+/epsilon1/prime i,∀i. (33) 15 In the matrix form, (33) can be rewritten as ξ=Hθ+/epsilon1, (34) where ξ,/bracketleftbigg ξ1ξ2... ξNT,/bracketrightbigg| ,θ,/bracketleftbigg δxn+1δyn+1/bracketrightbigg| ,/epsilon1,/bracketleftbigg /epsilon1/prime 1/epsilon1/prime 2... /epsilon1/prime NT/bracketrightbigg| and H,1 c ˆxn−x1 ˆRn1ˆyn−y1 ˆRn1 ˆxn−x2 ˆRn2ˆyn−y2 ˆRn2...... ˆxn−xNT ˆRnNTˆyn−yNT ˆRnNT . The best linear unbiased estimator of θis given by [ ?] ˆθ= (H|H)−1H|ξ, (35) withthecovariancematrix Cˆθ=σ2 iI(H|H)−1,where I2istheidentitymatrix,and ˆθisbydeﬁnition/bracketleftbigg ˆδxn+1ˆδyn+1/bracketrightbigg| . Based on the received signal strength estimates, the BLUE may be redeﬁned (improved) as ˆθ= (H|WH )−1H|Wξ, (36) withWa diagonal matrix W, ω10... 0 0ω2... 0 ............ 0 0... ωNT (37) whereωi,ˆαi/summationtextη i=1ˆαiifˆαi∈Qandωiis zero otherwise. The current estimate of the receiver position (at time instant n+ 1) is ˆxn+1= ˆxn+ˆδxn+1, (38) ˆyn+1= ˆyn+ˆδyn+1. (39) 16 IV. Simulations A. Simulation Preliminaries For simulations, we choose the following parameter values [4]: Responsivity Rp= 0.4mA/mW, Lambertian order /lscript= 1,height of the transmitter h= 2.5meters, area of the photo detector Ar= 1cm2and the spectral density of Gaussian noise is N0/2 = 1.336×10−22W/Hz [4]. If the bandwidth of the light pulse Bis chosen to be 100 MHz (same as the modulation bandwidth of the LED), then this translates to a noise power of σ2=N0 2×2B= 2.672×10−14W.Such a modulation bandwidth also translates approximately to a minimum pulse duration of 1×10−8seconds (if we ignore the sidelobes). Hence, Tp≥10−8seconds. The peak power (intensity) of the transmitted signal r(t)in (6) ranges from 1 through 10 W/m2[3]. We use a rectangular pulse for the purpose of positioning: r(t) =A· 1[0,Tp](t),where 1E(·)is the indicator function over a measurable set E. The sampling interval ∆is set equal to 1/(2B)(inverse of the Nyquist rate) seconds. For the rectangular pulse, the sampling frequency Fsis2BwhereBis the eﬀective bandwidth of the pulse. The bandwidth Bcan be set equal large enough to cover the ﬁrst two or three sidelobes of thesincfunction (which is the Fourier Transform of the rectangular pulse). Hence, for example, ifBis to include the ﬁrst three sidelobes, then B=4 Tpwhich implies Fs= 2B=8 Tp. This leads to ∆ =Tp 8. The performance criterion for positioning accuracy is the mean Euclidean error between the true position and the estimated position of the receiver. This is deﬁned as the mean Euclidean distance (or the/lscript2norm) between the actual coordinates (x,y)of the position and the estimated coordinates (ˆx,ˆy).The Euclidean error is denoted by ||e||2where e,/bracketleftbigg (x−ˆx) (y−ˆy)/bracketrightbigg| .The mean of the Euclidean error is computed in a Monte Carlo fashion by simulating the experiment many number of times. For the weighted least squares algorithm, we have set the number of simulations equal to 5000 for the purpose of computing the means. However, for the case of BLUE based positioning, the step size of the movement of the user is small, and at the same time, we want the user to traverse the entire room/premises in order to estimate the “average” positioning performance. Hence, the number of simulations for BLUE algorithm is set much higher at 50000. For the weighted least squares estimates, we use a real number genetic algorithm in order to ﬁnd the minimum of (21). The population size for the algorithm is set at 200, and number of generations 17 at 400. The crossover and mutation rates are set at 0.75 and 0.50, respectively. The reader may refer to [22] for more details about the genetic algorithms. In order to run a given simulation for weighted least squares based positioning, we choose the true coordinate (x,y)from within the region of the room/location randomly (uniformly). Then, we produce estimates (ˆx,ˆy)based on the observations. For BLUE, we start at a random location in the room, and the quantities δxandδyare sampled independently of each other from a uniform distribution:U[−ν,ν]. We refer to the quantity νas the “step size.” The value of νis set at 0.2 meters for most simulations on BLUE positioning. We conduct simulations for diﬀerent kinds of LED arrangements in a given location, and for diﬀerent numbers of receiver sensors which are arranged in a symmetric pattern about the center of the receiver. The total number of sensors Mis set equal to 2, 4, 6 and 8 for the purpose of performance comparison. B. Simulation results and performance evaluation (weighted least squares) 1) Weighted least squares In Figure 5 and Figure 6, we represent the normalized histograms of the error (diﬀerence between thexcoordinate of the true position and its estimate using the weighted least squares approach) for diﬀerent LED conﬁgurations and numbers of receiver sensors for the WLS estimate. These histograms correspond to the following parameter values: r(t) = 6W for Figure 5 and Figure 6(a) and r(t) =is3W for Figure 6(b), Tpis10−8seconds,Tis4×10−8seconds,Rpis 0.4, lis 1,his 2.5 m,Aris0.012m2,N0is2×1.336×10−22W/Hz,Bis100/TpHz,σ0is√N0BW, Fsis2BHz,∆is1/Fsseconds and cis3×108m/s. The center of the square shaped receiver is sampled from a two dimensional uniform distribution in the square region 1[−3,3](x)· 1[−3,3](y).The four sensors are present at each of the four corners, and the distance of each sensor from the center is 0.0707 meters. We used 5000 measurements for the creation of these histograms. We now consider Figure 7 which shows the LED arrangement of the ﬁrst kind, and its correspond- ing mean Euclidean error performance. The LED’s are represented by blue dots, and the receiver sensors by red dots. We use the same values of the parameters that we have laid out in the second paragraph of Section IV-B, except that, now, we evaluate the positioning performance for a range of transmitted power levels. Moreover, we also analyze the eﬀect of total number of receiver sensors on the performance. We varied the transmitted power from 2 W/m2through 10 W/m2in steps of 18 Fig. 5. Histogram of the measurement error along x-axis for the LED conﬁguration shown in Figure 7(a) for the synchronous system. The total number of receiver sensors Mis 4 and the peak signal power is 6 Watts per LED. (a) 0.3 0.2 0.1 0.0 0.1 0.2 0.3 error (meters)0.02.55.07.510.012.515.017.5probability (b) Fig. 6. Histogram of errors along x-axis. Figure (a) corresponds to LED conﬁguration shown in Figure 8(a) and M= 6and peak power of 6 Watts per LED. Figure (b) corresponds to LED conﬁguration shown in Figure 9(a) with M= 2and peak power of 3 Watts per LED. Both plots refer to the synchronous systems. 2 W/m2.The area of the square shaped room is 36 m2.The value of ηwas set equal to 3 (we used the three best sources out of 4). The mean Euclidean error performance is evaluated for both the synchronous and asynchronous systems. As seen in Figure 7, the mean Euclidean error decreases as the LED power levels are increased. Moreover, as expected, the error performance improves as we increase the number of sensors on the receiver. We were able to achieve the mean Euclidean error of approximately 0.02 meters with 8 sensors when the peak transmitted power for each of the four LED’s was 10 Watts. Figure 8 shows another layout of the LED’s for the room of similar dimensions and area as discussed previously, and the mean Euclidean error performance for this arrangement. The value of ηwas set 19 (a) LED arrangement (axes dimensions in meters) 2 3 4 5 6 7 8 9 10 peak signal power per LED (Watts)101 100mean Euclidean error (meters)2 sensors sync 4 sensors sync 8 sensors sync 2 sensors async 4 sensors async 8 sensors async (b) mean Euclidean error performance Fig. 7. LED arrangement of the ﬁrst type and its mean Euclidean error performance. (a) LED arrangement (axes dimensions in meters) 2 3 4 5 6 7 8 9 10 peak signal power per LED (Watts)102 101 mean Euclidean error (meters)2 sensors sync 4 sensors sync 8 sensors sync 2 sensors async 4 sensors async 8 sensors async (b) mean Euclidean error performance Fig. 8. LED arrangement of the second type and its mean Euclidean error performance. equal to 3. These performance curves are expected to be better than the system in Figure 7 since we now have an extra LED source in the middle. The improvement in performance is more pronounced when both the SNR and the number of sensors is small. We were able to localize the receiver to within a mean Euclidean error of 0.01 meters for M= 8and peak power of 10 Watts. Figure 9 represents the LED arrangement of the third type in which we have a total of 8 LED’s arranged symmetrically in the room. The value of ηwas set equal to 4. We compare the mean Euclidean error of this system with the one shown in Figure 7 (containing 4 LED’s), and ﬁnd out that this system outperforms the system with 4 LED’s even thought the average power of the two 20 3 2 1 0 1 2 33 2 1 0123 (a) LED arrangement (axes dimensions in meters) 1.0 1.5 2.0 2.5 3.0 3.5 4.0 4.5 5.0 peak signal power per LED (Watts)102 101 100mean Euclidean error (meters)2 sensors sync 4 sensors sync 8 sensors sync 2 sensors async 4 sensors async 8 sensors async (b) mean Euclidean error performance Fig. 9. LED arrangement of the third type and its mean Euclidean error performance. systems is the same. With 8 sensors and a peak power of 5 Watts per LED, we were able to achieve the mean Euclidean error of around 0.007 meters. Since the time synchronization error aﬀects the time-of-arrival based positioning, we consider the eﬀect of synchronization error in Figure 10. For a pulse of width Tpof 10 nanoseconds, we plotted the error curves as a function of the synchronization error2standard deviation that ranges from 0.2 to 1.0 nanoseconds. The synchronization error is assumed to be a zero-mean Gaussian random variable. Surprisingly, the error performance of the system with larger number of transmitters suﬀers more as the synchronization error becomes large. Moreover, as expected, the asynchronous systems are immune to synchronization errors. C. Simulation results and performance evaluation (BLUE) Figure 11 shows the mean Euclidean error performance for the two LED conﬁgurations for the step sizeνof 0.2 meters. The other parameters concerning the histograms are the same as discussed in the second paragraph of Section IV-B. Figure 12 represents the mean Euclidean error of the two systems with diﬀerent LED arrangements as a function of signal power. We note that the system with larger number of LED’s outperforms the other system even though both systems have the same average power. This conclusion is similar to the the corresponding results obtained in the case of weighted least squares based positioning. 2This error corresponds to the error in the transmitter and receiver clock time