Next Article in Journal
SMU Open-Source Platform for Synchronized Measurements
Previous Article in Journal
Study on TLS Point Cloud Registration Algorithm for Large-Scale Outdoor Weak Geometric Features
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

Pose Estimation for Visible Light Systems Using a Quadrature Angular Diversity Aperture Receiver

by
Shengqiang Shen
1,2,*,
Jose Miguel Menéndez Sánchez
3,4,
Shiyin Li
1,2 and
Heidi Steendam
3
1
School of Information and Control Engineering, China University of Mining and Technology, Xuzhou 221000, China
2
Engineering Research Center of Intelligent Control for Underground Space, Ministry of Education, China University of Mining and Technology, Xuzhou 221000, China
3
Department of Telecommunications and Information Processing, Ghent University/IMEC, Interuniversity Microelectronics Centre, 9000 Ghent, Belgium
4
Facultad de Ingeniería en Electricidad y Computación, ESPOL Polytechnic University, Campus Gustavo Galindo Km. 30.5 Vía Perimetral, Guayaquil 90112, Ecuador
*
Author to whom correspondence should be addressed.
Sensors 2022, 22(14), 5073; https://doi.org/10.3390/s22145073
Submission received: 31 March 2022 / Revised: 23 June 2022 / Accepted: 24 June 2022 / Published: 6 July 2022
(This article belongs to the Section Physical Sensors)

Abstract

:
The quadrature angular diversity aperture (QADA) receiver, consisting of a quadrant photodiode (QPD) and an aperture placed above the QPD, has been investigated for pose estimation for visible light systems. Current work on pose estimation for the QADA receiver uses classical camera sensor algorithms well known in computer vision. To this end, however, the light spot center first has to be obtained based on the RSS. However, this is less straightforward than for camera sensors, as in contrast to such sensors where the relationships are linear, the RSS output from the QADA is a non-linear function of the light spot position. When applying closed form solutions or iterative methods for cameras on a QADA, the non-linearity will degrade their performance. Furthermore, since in practice the aperture is not always perfectly aligned with the QPD, a procedure to calibrate the receiver is needed. Current work on calibration requires additional sophisticated equipment to measure the pose during calibration, which increases the difficulty of implementation. In this paper, we target the above problems for pose estimation and calibration of the QADA receiver. To this end, we first study the effect of the strategy of differencing and normalization on the probability density function (PDF), a commonly applied strategy for the QPD’s robustness against RSS variation, and it is shown that the applied strategy results in a complex PDF, which makes an effective and efficient estimation hard to achieve. Therefore, we derive an approximated PDF in a simple closed-form, based on which the calibration and the pose estimation algorithms using the least squares principle are proposed. The proposed calibration does not require any information about the pose of the receiver and is robust to variation of the received power and imperfect knowledge of the radiation pattern of the LED, making it easy to implement. We also derive the corresponding Cramér-Rao lower bound on the misalignment to benchmark the performance of the misalignment and to serve as an indicator to determine the required signal-to-noise ratio (SNR) or number of LEDs to obtain a desired accuracy. The calibration and pose estimation are evaluated by means of a Monte Carlo simulation. Computer simulations show that this theoretical bound is close to the RMSE of the proposed estimator and that the proposed pose estimator outperforms the PnP algorithm.

1. Introduction

White LEDs can be modulated up to several MHz. Therefore, they can be used for wireless communication and positioning. As white LEDs are becoming the primary source of light, visible light positioning (VLP) has the potential to achieve positioning with low cost, low power consumption, and long lifetime [1]. Due to the dependency of the RSS values on the pose (position and orientation) of the receiver, it is expected that imperfect knowledge of the pose may strongly affect the performance of visible light communication (VLC) systems [2,3,4,5,6,7]. Recently, a compact receiver based on a quadrant PD (QPD) combined with a single aperture, i.e., the quadrature angular diversity aperture receiver (QADA), was proposed for VLP [8]. The QPD is a segmented photodiode consisting of four matched PDs on a single chip separated by very thin gaps. Combined with an aperture, the resulting receiver achieves angular diversity with a reasonably compact size. Furthermore, as only one aperture needs to be placed above the QPD, the receiver will only have a limited number of construction parameters, making it more robust to construction errors.
The light from an LED reaches the QPD through the aperture only and introduces a light spot of the same size and shape as the aperture. As the position of this light spot is determined by the position of the LED and the pose of the receiver, the overlap area between the light spot and the four quadrants contains information about the pose of the receiver. The RSS values in the four constituent PDs of the QPD are proportional to this overlap area, indicating that a comparison of these RSS values allows us to estimate the angle-of-arrival (AOA) for the different LEDs as well as the position and/or orientation (pose) of the receiver. The current work has been initiated to achieve this goal. Experiments in [8] show that high angular resolution in estimated AOA can be achieved. The positioning performance of the QADA receiver is investigated in several works [8,9,10,11]. The paper [12] considers 3D position estimation with 1D orientation estimation, i.e., it is able to determine the azimuth angle for a receiver that is pointed upwards. The paper [13] that focuses more on the communication aspects considers the effect of a tilt angle, which restricts the orientation to two degrees of freedom. Recently, the work [14] studied the joint 3D position and 3D orientation estimation, namely pose estimation, based on a classic non-iterative perspective-n-point (PnP) algorithm, a well known algorithm for camera sensors in computer vision. In their work, to employ the PnP, first the light spot center is obtained based on the RSS. However, this is less straightforward than for camera sensors, as in such sensors the relationships are linear, whereas the RSS output from the QADA is a non-linear function of the light spot position. Moreover, the strategy of differencing and normalization is commonly applied to the RSS, claimed by most works to improve the robustness against RSS variation, which further increases the non-linearity. These effects will degrade the performance of the PnP algorithm when applied to the QADA receiver; the PnP algorithm is sensitive to noise and as a result it lacks precision [15].
Furthermore, since in practice the aperture is not always perfectly aligned with the QPD, a procedure to calibrate the receiver is needed. The current work starts to consider the presence of misalignment between the aperture and the QPD. In [8], the authors manually adjusted the position of the aperture as a pre-step in the experiment, and the need for calibration is further emphasized in [9]. However, such a hardware calibration is a burdensome, time-consuming procedure. A better approach is to analyze the effect of the aperture misalignment and compensate for it through signal processing. For example, in [12], the authors propose a calibration method to compensate for the mismatch using signal processing. As the authors assume in this calibration phase that the relative distance between the LED and the receiver is known, additional equipment such as Optitrack is required to determine the ground truth on this distance. This limits the usability of the considered method.
In this paper, we target the above problems for pose estimation and calibration of the QADA receiver. On the one hand, we propose a calibration procedure using signal processing that does not require the knowledge of the pose of the receiver; thus no additional hardware is needed to obtain the ground truth on this pose. On the other hand, a pose estimation algorithm is proposed that directly estimates the pose from the (normalized) RSS values, achieving better performance compared to the method that estimates the light spot center and then applies the camera algorithm. To this end, we first investigate the strategy of differencing and normalization to the RSS, i.e., where we only look at the relative differences between the normalized RSS values to obtain an observation that is robust against variation in received power and offsets common to all quadrants [12]. We investigate the distribution of the observation obtained from the noisy RSS values, based on which we then propose the algorithm to estimate the pose of the receiver that is more robust to noise.
The main contributions of this paper are as follows:
  • We first model the RSS vector of the QADA at an arbitrary pose using the perspective projection model. With the help of this model, we derive an explicit expression relating the RSS with the intrinsic parameters, i.e., the aperture height and the misalignment of the aperture, and the extrinsic parameters, i.e., the position and orientation of the receiver.
  • We use the strategy of the normalized differences of the RSS values to make the estimator robust against variations in the transmitted power and radiation pattern. To be able to derive a simple estimation algorithm, we replace the complex true PDF of the resulting observations by an approximated Gaussian PDF based on the first order Taylor series approximation of the observation.
  • Using this approximated PDF, we propose a calibration algorithm based on the least squares (LS) principle. The algorithm jointly estimates the intrinsic and extrinsic parameters from which the intrinsic parameters are extracted. This estimation is performed in an iterative way, where the principle of optimization on manifolds is used. After calibration, a simplified version of the estimation algorithm is used to extract the pose, as now we use the calibrated intrinsic parameters as prior knowledge for the misalignment.
  • To evaluate the optimality of the proposed algorithms, we derive the misspecified Cramér-Rao bound (MCRB) to take into account the effect of the approximated PDF. We compare the MCRB with the Cramér-Rao bound (CRB) for the detected RSS values to quantify the performance loss due to using the normalized differences of the RSS values to make the estimator robust against imperfect knowledge of the transmitted power and radiation pattern. The designed algorithms are verified by Monte Carlo simulations.
The rest of the paper is organized as follows. The channel link model and the expressions for the observation and output vectors are presented in Section 2. The approximated PDF and the calibration system model are provided in this section as well. The calibration and the pose estimation algorithms using the LS principle are proposed in Section 3. Subsequently, the theoretical bound is derived in Section 4. The Monte Carlo numerical comparison is given in Section 5. Finally, some concluding remarks are given in Appendix A.
Notation: Scalars are denoted in italic, e.g., x. Lowercase boldface indicates a column vector, e.g., x . Uppercase boldface denotes a matrix or a set, e.g., X , with I N representing an N × N identity matrix, and 0 N × M representing an N × M zero matrix. The vector e i is a unit basis vector with its i th element being 1, and the operator · ¯ converts the Cartesian coordinates x into the homogeneous coordinates x ¯ = [ x T , 1 ] T . Matrix transpose and inverse are indicated by superscript T , and 1 , respectively. The Euclidean norm is denoted by · , and the expectation operator is denoted by E { · } . x = / x denotes the Del operator, while Δ x z = z x T denotes the Hessian operator. The rectangular function Π · is defined as
Π x = Δ 1 , x 1 . 0 , x > 1 .
The set of all real numbers is denoted by R . The group of all rotation matrices, i.e., the special orthogonal group, is denoted by S O ( 3 ) , and the associated Lie algebra is denoted by so ( 3 ) . The special Euclidean group is denoted by S E 3 . The operator · converts a 6 × 1 vector into a member of se 3 – the Lie algebra of S E ( 3 ) – which is,
a T , b T T = b × a 0 T 0 R 4 × 4 , a , b R 3 × 1 ,
while the operator · × converts a 3 × 1 vector into so 3 , which is,
a 1 a 2 a 3 × = 0 a 3 a 2 a 3 0 a 1 a 2 a 1 0 .
The operators · and · are, respectively, the inverse operators of · and · × . Additionally, exp · and log · are the matrix exponential and matrix logarithm functions, respectively.

2. System Description

In this section, we describe a visible light system that is able to detect signals transmitted by different LEDs from which the pose of the receiver can be determined. To avoid interference between the reference signals for pose estimation, we assume the system adopts multicarrier multiplexing based on (DC-biased) sinusoidal waves [16,17], where a portion of the subcarriers is assigned to the LED anchors to be used as reference signals for pose estimation, while the rest of the subcarriers can be used for communication. Each anchor is assigned a different subcarrier, implying that due to the orthogonality of the subcarriers, the reference signals from the different LEDs can be separated in a simple way. We assume the receiver has knowledge of the used subcarrier and position of each LED, e.g., because this information is contained in the reference signal or is available in a database accessible by the receiver.
The rest of this section is organized as follows. We first model the channel link between a single LED and the QADA and give the expressions for the output vector and observation. We further derive the PDF of the observation. As this PDF is too complex for practical estimation purposes, we also derive a simple closed-form approximation. Finally, we design a system model for calibration and discuss the goals of the pose estimation problem.

2.1. Receiver Structure

The considered visible light system uses a QPD consisting of four matched photodiodes deposited on a single chip separated by very thin gaps. To convert the incident light of an LED into a light spot on the QPD, we place an opaque screen with a circular aperture at a specific height above and parallel to the surface of the QPD [8] as shown in Figure 1 (a square aperture has the disadvantage of being sensitive to rotational offsets, which means that the rotational offset must be calibrated and compensated for, while a circular aperture is circularly symmetric, implying it is inherently robust to rotational offsets. Therefore, in this paper, we focus on the circular aperture). The light spot on the QPD’s surface, as shown in Figure 2, leads the QPD to produce a unique quad-tuple of RSS values proportional to the overlap areas, i.e., the RSS output is uniquely determined by the position of the light spot. Because of this, the position of the light spot can be estimated from the RSS values, which in turn, e.g., is useful in detecting an object’s angle coordinates [8].
The aperture is chosen so that its diameter is very large compared to the wavelength of the light. Further, as the distance between the LED and the QPD is very large compared to the wavelength, the receiver is in the far-field region of the LED, implying the incident waves can be considered as plane waves. As a result [18], the only light that reaches the QPD is the light that passes through the aperture, and the incident light will introduce a circular light spot of the same size as the aperture on the plane of the QPD. The position of the light spot in the QPD plane not only depends on the position of the light source, but also on the pose of the receiver, implying the RSS outputs in the constituent PDs of the QPD are functions of the position of the light source and the pose of the receiver. We assume that the light spot overlaps with all quadrants and completely falls within the active area of the QPD (in case the light spot does not overlap with all quadrants, one or more of the RSS values will be zero. This will obstruct the strategy of differencing and normalization of the RSS values, which is used to make the estimation algorithm robust against variations in the transmitted optical power. In case the light spot partially falls off the QPD, the relationship between the RSS values and the pose of the receiver becomes very complex, implying the resulting estimation algorithm will have high complexity). To ensure that the first condition is fulfilled, we assume the receiver considers only LEDs that lead all quadrants to produce a sufficiently high RSS output. The second condition implies that the aperture diameter is upper limited by (roughly) the dimension of a quadrant, i.e., half the diameter of a circular QPD or half the side of a square QPD. This guarantees a simple unique relationship between the quad-tuple of RSS values generated by the QPD and the pose of the receiver, and by combining the information from different light sources satisfying the above conditions, the pose of the receiver can easily be determined.
In the following, we will describe the relationship between the pose of the QADA and the position of the light spot, and then we will derive the relationship between the light spot and the RSS output.

2.2. Channel Link Model

To obtain the quad-tuple of RSS values as a function of the pose, we note that the dependency of the RSS values on the pose is twofold. Firstly, the channel gain is determined by the radiation pattern of the LED and the irradiance attenuation, implying the incident power at the QPD (and thus the optical power in the light spot) depends on the relative position of the light spot with respect to the LED. Secondly, the RSS values in the quadrants are proportional to the overlap areas of the light spot with the quadrants of the QPD, which is determined by the position of this light spot with respect to the center of the QPD. To express the RSS values as function of the pose of the receiver, we define the following three coordinate systems: (1) The first coordinate system, i.e., the system frame, corresponds to the inertia frame and is used to describe the pose of the receiver with respect to the LEDs. Furthermore, we define two coordinate systems that are attached to the receiver and will be used to describe the position of the light spot and misalignment due to construction errors; (2) In the first receiver attached coordinate system, i.e., the receiver frame, the x-y plane is the plane of the aperture, and the z-axis is the normal to that plane. The origin O r of this receiver frame coincides with the center of the aperture. We will use an extra subscript r to indicate coordinates represented in the receiver frame, i.e., the coordinates x r in the receiver frame correspond to their counterpart x in the system frame. This coordinate system is used to find the distance vector of the light spot to the LED as a function of the pose; (3) The second receiver attached coordinate system is the QPD frame, where the x-y plane is the surface of the QPD, the origin O q is placed in the center of the QPD, and the x and y axes are aligned with the gaps between the quadrants. We will use an extra subscript q to indicate coordinates represented in the QPD frame, i.e., the coordinates x q in the QPD frame correspond to their counterpart x in the system frame. Although this coordinate system is closely related to the receiver frame, its benefit lies in the simplification of the computations, as it enables us to convert the 3D receiver frame coordinates to the 2D coordinates in the QPD frame (we assume that the aperture plane is parallel to the plane of the QPD, and the x- and y-axes of the QPD and receiver frames are aligned. This implies that the transformation of vectors in the aperture frame into vectors in the QPD frame is a simple translation). This approach is commonly used in the perspective projection model, which we employ to determine the position of the light spot in the QPD plane, and from which we compute the overlap areas of the light spot and the QPD. In the following, we will determine the coordinates of the light spot generated by an LED.
Assume that in the system frame, the pose of the receiver is described by its orientation R S O ( 3 ) , i.e., R belongs to the Special Orthogonal group and position r R 3 × 1 , where the position of the receiver is defined as the position of the center of the aperture. We further assume the LED has position r L and normal vector n L in the system frame. To calculate the position of the light spot in the QPD frame, we notice that the line going through the position r L of the LED and the centers r of the aperture and r S R 3 × 1 of the light spot can be seen as a ray in the perspective projection model of a center of projection located at the center of the aperture and having focal length equal to the aperture height | h | , i.e., the distance between the planes of the aperture and the QPD (see Figure 3). The perspective projection model allows us to relate the 3D coordinates of the LED—the `object’—in the system frame to the 2D coordinates of the light spot —the `image’—in the QPD frame. Note that the position of the light spot in the QPD frame is specified by the ( x , y ) coordinates only. To make the paper self-contained, we include the background on the perspective projection model. Taking this into account and using the projection principle, the position r S , q R 2 × 1 of the center of the light spot in the QPD frame is given by [19]:
r S , q = k m , q + e 1 T e 2 T I 3 , 0 h e 3 T T r ¯ L T r ¯ L A B C .
The matrix T S E 3 is defined through its inverse
T 1 = R r 0 1 × 3 1 .
Furthermore, the misalignment vector k m , q expresses the errors in the alignment of the aperture with respect to the QPD. Due to construction errors, a small horizontal offset between the origins of the QPD frame and the receiver frame may occur. This is illustrated in Figure 3, where the brown dot does not coincide with the center O q of the QPD. The horizontal misalignment vector, which takes into account this offset in the transformation of the light spot coordinates in the receiver frame to the coordinates in the QPD frame, is defined as k m , q [ u m , v m ] T , with u m and v m the misalignment in the x and y direction, respectively. Taking into account (4), we see that the projection principle consists of three steps:
  • Step A: This step converts the LED coordinates r ¯ L in the system frame to the coordinates r ¯ L , r = T r ¯ L in the receiver frame using the matrix T .
  • Step B: This scales r ¯ L , r along the projection line with the factor h h L where h L = e 3 T T r ¯ L and gets rid of the last component of the homogeneous coordinates through the 3 × 4 matrix I 3 , 0 , consisting of the 3 × 3 identity matrix extended with a zero column, as this last component is irrelevant for the determination of the position of the light spot. This results in the position r S , r = I 3 , 0 h e 3 T T r ¯ L r ¯ L , r of the light spot in the receiver frame.
  • Step C: This transforms r S , r into the 2D coordinates r S , q by discarding the z coordinate through the mapping of r S , r on the 2 × 3 matrix e 1 T e 2 T and adds the misalignment vector k m , q .
Equation (4) can be rewritten in the following compact form:
r S , q = K T r ¯ L e 3 T T r ¯ L .
The matrix K , which is used in the perspective projection model, is called the intrinsic parameter matrix and describes the internal properties of a receiver, i.e., that are independent of the pose. It only depends on the misalignment vector k m , q and the aperture height h:
K = h 0 u m 0 0 h v m 0 .
On the other hand, the matrix T contains the pose information of the receiver and is independent of the internal properties of the receiver. Therefore, it is called the extrinsic parameter matrix. Similarly, taking into account (4), we can write the position r S , r of the center of the light spot in the receiver frame (in homogeneous coordinates) in the following compact form:
r ¯ S , r = K 0 T r ¯ L e 3 T T r ¯ L with K 0 = h I 3 0 e 3 T 0 ,
and the position r ¯ S of the light spot in the system frame (in homogeneous coordinates) is obtained through:
r ¯ S = T 1 r ¯ S , r .
To determine the RSS values in the quadrants of the QPD, we first need to compute the overlap area between the light spot and the quadrants of the QPD. Let us assume the aperture has radius l. Taking into account we already assumed that the light spot completely overlaps with the QPD, i.e., no part of the light spot falls outside the QPD, and that all quadrants overlap with the light spot, the overlap areas can easily be computed based on the position r S , q of the center of the light spot in the frame of the QPD. Defining the overlap area A i > 0 of the ith quadrant, i = 1 , , 4 and the angles α j [ 0 , π ] with [ r S , q ] j = l cos α j , j = 1 , 2 , and the operator [ r S , q ] j as the jth element of r S , q (see Figure 2), we obtain the overlap areas by solving the following set of equations:
  • The light spot has area π l 2 , implying A 1 + A 2 + A 3 + A 4 = π l 2 .
  • Taking into account that the area of a circular segment with central angle 2 α , α [ 0 , π ] equals l 2 2 ( 2 α sin 2 α ) , we find that A 2 + A 3 = ( α 1 sin α 1 cos α 1 ) l 2 and A 3 + A 4 = ( α 2 sin α 2 cos α 2 ) l 2 .
  • The overlap area A 1 with the first quadrant is a combination of a circle sector with central angle 3 π 2 ( α 1 + α 2 ) having an area 3 π 2 ( α 1 + α 2 ) l 2 2 , a rectangle with area [ r S , q ] 1 [ r S , q ] 2 = l 2 cos α 1 cos α 2 and two triangles with areas l 2 2 cos α j sin α j , j { 1 , 2 } .
Note that, depending whether cos α j > 0 (if α j [ 0 , π 2 ] ) or cos α j < 0 (if α j [ π 2 , π ] ), the above contributions will combine positively or negatively, resulting in the wanted overlap area A i . Hence, the overlap area vector A = [ A 1 , , A 4 ] T is the solution of M a A = s with
M a = 1 0 0 0 0 1 1 0 0 0 1 1 1 1 1 1 and s = l 2 3 π 4 α 1 + α 2 2 + sin 2 α 1 + sin 2 α 2 4 + cos α 1 cos α 2 α 1 sin α 1 cos α 1 α 2 sin α 2 cos α 2 π .
Now that we obtained the overlap areas of the light spot with the PDs of the QPD and related them to the pose of the receiver, the RSS values in the different PDs can be determined. To this end, we define the distance vector between the transmitter and the light spot as v = r S r L . With this definition, we can write the distance v between the LED and the light spot, and the incidence angle θ , i.e., the angle between the vector v and the normal n Q = Re 3 of the QPD as:
v = v ,
cos θ = T 1 e 3 T v ¯ v ,
where the last equation holds due to the equalities n Q T , 0 T = T 1 e 3 , obtained from (5), and n Q T v = n Q T , 0 v ¯ . The channel gain h i for the i th quadrant of the QPD is given by [20]
h i = Γ ( n L , v ) cos θ Π θ 2 π v 2 A i ,
where Γ ( n L , v ) is the radiation pattern of the LED pointing n L and is evaluated at direction v , Π θ = Π ϕ / ϕ F O V . The factor Π θ in (13) implies that the QPD can detect the light only when the LED is within its FOV, i.e., 0 ϕ ϕ F O V . Taking into account (11)–(13), the channel gain vector is given by
h k , T = Π θ Γ n L , v T 1 e 3 T v ¯ 2 π v 3 A ,
and is a function of the intrinsic parameter vector k [ u m , v m , h ] T and the pose T . From this channel gain, we obtain the RSS output vector y = [ y 1 , , y 4 ] T in the four PDs:
y = g + η ,
where g = R p P t h , with R p the responsivity of the QPD, and P t the transmitted power. Due to coupling between the quadrants, the noise components are correlated Gaussian distributed random variables η N ( 0 , σ η 2 C ) , where C is the correlation matrix with C i , i = 1 and C i , j = ρ if i j , ρ is the correlation between the noise of different quadrants, which can be determined by comparing the noise values in the different quadrants during an experiment and calculate the resulting correlation, and σ η 2 is the noise variance.
In the derivation of g in (15), it is assumed that the transmitted optical power is perfectly known, and the LED is a perfect Lambertian radiator. However, directly estimating the intrinsic and/or extrinsic parameters from the RSS output vector will result in inaccurate estimates. To make the estimator more robust to changes in the incident power, the normalized differences between the RSS values along the two axes can be used, as
t x = ( y 1 + y 4 ) ( y 2 + y 3 ) y 1 + y 2 + y 3 + y 4 t y = ( y 1 + y 2 ) ( y 3 + y 4 ) y 1 + y 2 + y 3 + y 4 .
The approach (16) will make the observations insensitive to variations in the transmitted optical power and radiation pattern provided that the RSS values are noise-free or affected by a common offset so that the normalized differences of (16) are noiseless observations. However, the noise induced by ambient light in visible light systems is non-negligible, and due to the randomness of the shot and thermal noise, the noise present in the different quadrants is in general not equal. The presence of this non-identical random noise will have an impact on the distribution of t x and t y . As the estimation of the intrinsic and extrinsic parameters relies on proper knowledge of this distribution, we need to evaluate the effect of the noise on the distribution to see if the observations (16) are still insensitive to variations in the transmitted power when noise is present. Therefore, in the following section, we will analyze the distribution of the normalized differences t x and t y .

2.3. Approximation to the PDF of the Observation

To find the distribution of the normalized differences, we first rewrite the numerator and denominator of t x and t y in a more compact form by introducing m m 1 , m 2 , m 3 T = M y g and w w 1 , w 2 , w 3 T = M y η with
M y = 1 1 1 1 1 1 1 1 1 1 1 1 ,
implying (16) can be rewritten as t = m 1 + w 1 m 3 + w 3 , m 2 + w 2 m 3 + w 3 T . Taking into account that the components of w are linear combinations of Gaussian distributed random variables, w has a multivariate normal distribution, i.e., w N 0 , Σ w with
Σ w = 4 σ η 2 1 ρ 0 0 0 1 ρ 0 0 0 1 + 3 ρ .
Note that the numerator m 3 + w 3 = i y i corresponds to the sum of the RSS values in the different quadrants. Hence, m 3 and σ w 3 2 4 1 + 3 ρ σ η 2 are the average and variance of the total received signal strength, which implies we can define ξ m 3 2 / σ w 3 2 as the received signal-to-noise ratio (SNR) at the receiver. Furthermore, we define
μ m 1 m 3 , m 2 m 3 T ,
i.e., the ratio of the averages of the denominator and numerator of t . Obviously, in the absence of noise, the observation equals t = μ . We will show that:
  • μ is independent of the channel gain parameters and the transmitted power.
  • At high SNR, the distribution of t is given by t ˙ N ( μ , Σ ) with
    Σ = ξ 1 1 ρ 1 + 3 ρ I 2 + μ · μ T .
Let us first take a closer look at the dependency of μ on the channel gain and the transmitted power. Defining h ( k , T ) = λ A , and taking into account that the factor λ is independent of the considered quadrant, it follows that μ is independent of the channel gain parameters and the transmitted power contained in the factor λ , but only depends on the overlap area vector A . Substituting the solution of M a A = s in (19), we find after some straightforward computations that
μ i = 1 2 α i cos α i sin α i π .
Hence, the vector μ solely depends on the angles α i that are determined by the coordinates of the center of the light spot in the QPD frame. These coordinates, and therefore also μ , only depend on the intrinsic and extrinsic parameters through (6) but not on the channel gain parameters nor the transmitted power.
To obtain the distribution of t at high SNR, we first notice that t x = m 1 + w 1 m 3 + w 3 and t y = m 2 + w 2 m 3 + w 3 are both a ratio of two independent normal variables with non-zero means, and as they share the same denominator, they are correlated. The distribution of a single ratio of normal variables has been investigated in the literature, e.g., [21,22,23]. In these works, it is shown that a closed-form expression for the PDF in the general case for such a ratio is very complex as its shape differs significantly with its parameters, i.e., in some cases it resembles the Cauchy distribution, in other cases a normal distribution or a bimodal distribution. Several closed-form approximations were discussed in the literature, e.g., the approximate normal distribution. In this paper, we need the distribution of two correlated ratios of normal variables. As this is an even more complicated situation compared to the single ratio case, it is clear that the closed-form expression for the PDF will be even more complex. To simplify the analysis, we first notice that at high SNR, w i m 3 1 . Expanding t with respect to w using the Taylor series, keeping up to the linear terms in w , we obtain:
t μ + M w w with M w = 1 m 3 0 m 1 m 3 2 0 1 m 3 m 2 m 3 2 .
This first-order approximation directly leads to the distribution t ˙ N ( μ , Σ ) with Σ = M w Σ w M w T . Substituting M w and Σ w (18) in Σ , and using the definition of the SNR ξ , it follows that the covariance matrix Σ reduces to (20).
Following the above analysis, at high SNR we can design an estimator based on μ to achieve robustness against variations in the transmitted optical power and radiation pattern to estimate parameters from t .

3. Calibration and Pose Estimation

Using the approximate PDF of the normalized differences derived in the Section 2, we can now estimate the pose of the receiver, taking into account the relations (19) and (21). When the intrinsic parameters are known, the estimation of the pose, i.e., the extrinsic parameters, is straightforward. However, in practice, the intrinsic parameters are not prior known and need to be estimated because inaccurate knowledge of the misalignment parameters may result in significant biases in the pose estimates. As these intrinsic parameters can be considered fixed once the receiver is assembled, the intrinsic parameters can be determined once during a calibration process. The standard approach is to determine the intrinsic parameters assuming the pose of the receiver is known. However, this approach requires accurate knowledge of the pose of the receiver. Although it is possible to accurately determine the pose of the receiver, this requires costly equipment, in particular to determine the orientation of the receiver, and a laborious procedure to determine the pose of the receiver before each measurement. Therefore, we propose in this section a method to estimate the intrinsic parameters not requiring the knowledge of the pose, making it simpler and less costly than the standard calibration procedure. Note that our approach differs from the classical algorithm commonly used in camera-based systems, e.g., the non-iterative PnP algorithm and its iterative version [24], as in our method we directly estimate the pose from the normalized differences, while in the classical algorithm used for cameras, first the position of the light spot is determined, after which the pose is extracted. Although the algorithm used for cameras can be adapted to solve the problem at hand, we expect that this algorithm will perform worse than the algorithm proposed in this paper, as compared to the high resolution pixel measurement from the camera sensor, the RSS output from the QADA typically has a lower SNR. Furthermore, because of the non-linear relationship between the center of the light spot and the RSS values, the algorithm will be suboptimal for the considered system [15,25]. As the proposed method directly estimates the pose from the RSS values, it also differs from classical techniques used to determine the position, i.e., where first the AOA or distances to the LEDs are estimated before determining the position using, e.g., trilateration or triangulation.

3.1. Calibration Procedure

The calibration procedure proposed in this paper makes use of the setup shown in Figure 4, where N L LEDs are installed on a plane in the system frame x O s y , and the positions of these LEDs are fixed and known. The receiver observes these LEDs at N T different poses. Hence, the receiver observes for each LED N T quad-tuples of RSS values. The resulting observations are arranged in the 4 N T N L × 1 vector y ˘ R 4 N T N L × 1 :
y ˘ = g ˘ + η ˘ ,
where the vector η ˘ N 0 , Σ ˘ η collects the noise, and Σ ˘ η = σ η 2 I N T N L C (see Section 2.2 for the definition of the noise correlation matrix C ). The vector g ˘ can be rewritten as g ˘ = R p P t , 1 g 1 , 1 T , P t , 2 g 1 , 2 T , , P t , N L g N T , N L T T , where g i , j R 4 × 1 is the channel gain between the j th LED and the receiver at the i th pose T i 1 , and P t , j is the transmitted power of the j th LED. For each LED and each pose, we compute the normalized differences and arrange them in the 2 N T N L × 1 vector t ˘ R 2 N T N L × 1 . Taking into account that the noise in the different observations is statistically independent, we can write the Gaussian approximation of the PDF of t ˘ as:
p m t ˘ = 1 det ( 2 π Σ ˘ ) e ( 1 2 t ˘ μ ˘ Σ ˘ 2 )
where μ ˘ = μ 1 , 1 T , μ 1 , 2 T , , μ N T , N L T T . The components μ i , j correspond to the j th LED and the i th pose and are defined through (19), while the covariance matrix Σ ˘ is the block diagonal matrix given by
Σ ˘ = diag ( Σ 1 , 1 , Σ 1 , 2 , , Σ 1 , N L , Σ 2 , 1 , , Σ N T , N L )
In this covariance matrix, Σ i , j is the covariance matrix for the noise of t i , j (see (20)). These covariance matrices Σ i , j are a function of μ i , j , which in turn depend on the parameters to be estimated. As a consequence, maximum likelihood estimation will be complex. Therefore, we will use the least squares (LS) method to estimate the intrinsic and extrinsic parameters. Let us define the parameter set Θ = { k , T 1 , , T N T } , then the LS estimate is given by:
Θ ^ = arg min Θ 1 2 t ˘ μ ˘ 2 , s . t . R i T R i = I , det R i = + 1 .
The constraints R i T R i = I and det R i = + 1 imply that R i , which defines the orientation of the receiver at pose i and is enclosed in the transform matrix T i 1 S E 3 (see (5)), is a rotation matrix that is a member of S O 3 . Unfortunately, (26) has no closed-form solution, implying we need to resort to an iterative algorithm to estimate Θ . As the solution of the above (iterative) constrained optimization problem is complex and cumbersome in Euclidean space, we first notice that S E 3 is a manifold and convert the above optimization problem into an unconstrained optimization problem on manifolds [26]. In what follows, we use the Gauss–Newton method on manifold S E 3 to solve (26), similarly as in [27]. At each iteration, the update direction is calculated by
Δ k , Δ T 1 , , Δ T N T T = Θ μ ˘ t ˘ μ ˘
where ( · ) is the Moore–Penrose pseudoinverse, and Θ μ ˘ = [ Θ μ 1 , 1 , Θ μ 1 , 2 , , Θ μ N T , N L ] T R 2 N T N L × 6 N T + 3 the gradient of μ ˘ with respect to Θ . The derivation of the components Θ μ i , j can be performed in a similar way as in [27], and the result is given in Appendix A. As the update directions (27) depend on k and all T i , the intrinsic parameters and the poses need to be estimated jointly:
k t + 1 = k t + τ k Δ k ( k t , T 1 t , , T N T t ) T i t + 1 = exp τ T Δ T i ( k t , T 1 t , , T N T t ) T i t , i { 1 , , N T }
where τ k and τ T control the incremental step size for k and T i , respectively.
The Gauss–Newton method needs an initialization Θ 0 to start the iterative estimation process. To obtain this initial estimate, we consider the direct linear transformation method [24]. Assuming the positions of at least four LEDs, and the positions of the light spots for these LEDs are known, the direct linear transformation method gives a closed-form coarse estimate for Θ 0 . Although in our problem, the positions of the light spots are not known, we can estimate them from the normalized differences t i , j , which is a noisy version of μ i , j (22). This μ i , j is a function of the angles α , = 1 , 2 (21), which in turn are related to the position r S , q , i , j of the light spot, i.e., we can write μ i , j as a function of r S , q , i , j : μ i , j ( r S , q , i , j ) . In many situations, the noise is relatively small, implying we can neglect the presence of the noise, i.e., t i , j μ i , j . Inverting the function μ i , j ( r S , q , i , j ) , we obtain the coarse estimate r ^ S , q , i , j R 2 × 1 of the light spot for the j th LED and the i th pose:
r ^ S , q , i , j = μ i , j ( 1 ) ( t i , j ) ,
where ( · ) ( 1 ) is the inversion operator. Unfortunately, due to the non-linearity of μ as a function of its argument r S , q , no closed-form expression can be found for this inverse. However, as μ monotonically increases with r S , q , we can precompute the function μ for a set of values r S , q and save them in a look-up table. Based on this table, the inverse can be determined by interpolating between the values available in the look-up table.

3.2. Theoretical Lower Bound

In the Section 3.1, we considered the LS estimation of the intrinsic parameters based on the approximated PDF of t ˘ in the presence of the (unwanted) extrinsic parameters T i . To gain insight into the optimality of the designed calibration algorithm and to investigate how the SNR, N L or N T would effect the accuracy of the estimator, in this section we derive the theoretical lower bound for the mean squared errors of the intrinsic parameters h and k m , q . Since the estimation includes the unwanted parameter T i and since the estimator is designed based on an approximated PDF, we will derive the misspecified Cramér-Rao bound (MCRB) [28] for the whole parameter set Θ and only keep the left upper 3 × 3 submatrix corresponding to the intrinsic parameters.
The MCRB for the whole parameter set Θ is given by
MCRB ( Θ ) = ( M 1 ( Θ ˚ ) ) 1 M 2 ( Θ ˚ ) ( M 1 ( Θ ˚ ) ) 1 + Bias ( Θ ˚ , Θ ) .
In (30), Θ ˚ is the parameter set to which the estimation converges due to the PDF approximation assumed by the estimator, given by
Θ ˚ = arg min Θ D p ( t ˘ ) | | p m ( t ˘ | Θ ) ,
where D p ( t ˘ ) | | p m ( t ˘ | Θ ) is the Kullback–Leibler divergence (KLD) between the true ( p ( t ˘ ) ) and the approximated ( p m ( t ˘ | Θ ) ) PDFs. The matrices M 1 , M 2 , and Bias ( Θ ˚ , Θ ) are, respectively, given by [28]
M 1 = E t ˘ Δ Θ Θ ln p m ( t ˘ | Θ ) ,
M 2 = E t ˘ Θ T ln p m ( t ˘ | Θ ) Θ ln p m ( t ˘ | Θ ) ,
Bias ( Θ ˚ , Θ ) = ϵ ( Θ ˚ , Θ ) ϵ ( Θ ˚ , Θ ) T ,
where ϵ ( Θ ˚ , Θ ) is the error vector between Θ ˚ and Θ . The error vector is defined as ϵ ( Θ ˚ , Θ ) = ϵ k T , ϵ T 1 T , , ϵ T N T T T with ϵ k = k ˚ k and ϵ T i = log ( T i T ˚ i 1 ) . After a few rearrangements, it can be seen that the expectations in (32) and (33) depend on p ( t ˘ ) through its mean and covariance only. Due to the complexity of the true PDF p ( t ˘ ) , this mean and average must be calculated numerically via Monte Carlo integration.
To analyze the effect of N T , N L and SNR, we take in account that since the approximated PDF p m ( t ˘ | Θ ) is Gaussian, the minimum of the KLD (31) is obtained when the first two moments of the true and the approximated PDFs are matched [29]. When the SNR is high, the first two moments of the approximated PDF approach those of the true PDF. As a result, the estimator is asymptotically unbiased, i.e., Bias ( Θ ˚ , Θ ) 0 for SNR 1 , so that M 1 M 2 , MCRB ( Θ ) M 2 1 and
M 2 i N T j N L Θ T μ i , j Σ i , j 1 Θ μ i , j .
It can be seen that when N T or N L increases, the number of terms in (35) increases. As the matrices being summed in (35) are symmetric and positive-semidefinite, this implies that M 2 is non-decreasing in the sense of the Loewner order. Furthermore, as Σ i , j 1 is proportional to the SNR, increasing the SNR also enlarges the Loewner order of these matrices being summed and thus leads the partial order of M 2 to enhance. As a result, the MCRB is a non-increasing function of N T , N L , or SNR. While an increase of N T , N L , or SNR improves the robustness of the estimator against noise, these parameters do not have an impact on the feasibility of the estimation. The number of observations required to make proper estimation feasible follows from the analysis from [24] dealing with intrinsic parameter estimation for cameras. Translating the results of that paper to the problem at hand, it follows that observations at two different poses provide sufficient constraints to avoid depth ambiguity, and at each pose, at least four LEDs should be observed to have enough information for the calibration.
In this paper, we estimated the intrinsic parameters from the normalized differences instead of the quad-tuple of RSS values, where the PDF of the normalized differences was approximated by a Gaussian distribution. To evaluate the effect of both the approximation and normalized differencing, we compare the misspecified Cramér-Rao bound MCRB Θ (30) with the Cramér-Rao bound for the quad-tuple of RSS values y ˘ observed by the QPD, i.e., CRB Θ , corresponding to the true PDF p ( y ˘ ) . As y ˘ is Gaussian distributed, the derivation of the Fisher information matrix J Θ and thus CRB ( Θ ) = J 1 ( Θ ) [30] is straightforward and yields
J Θ = Θ g ˘ T Σ ˘ η 1 Θ g ˘ ,
where Θ g ˘ = Θ g 1 , 1 , Θ g 1 , 2 , , Θ g N T , N L T R 4 N T N L × 6 N T + 3 denotes the gradient of g ˘ with respect to Θ .

3.3. Pose Estimation

Once the intrinsic parameters are estimated during the calibration process, the estimated parameters can be used for accurate pose estimation of the receiver. To estimate the pose, a similar procedure can be used for the calibration by means of the Gauss–Newton method, except that now it is assumed that the intrinsic parameters k are known. At pose i, the observation t i only depends on the pose T i and is independent of other poses. Therefore, we can estimate the poses independently, in contrast to the calibration procedure that needed to consider the observations of multiple poses due to the presence of the unknown intrinsic parameters. Again, the estimation is obtained with an iterative procedure of which the update step is the same as in (28):
T i t + 1 = exp τ T Δ T i ( k , T i t ) T i t
Similar to the calibration procedure, we can obtain a coarse estimate by means of the direct linear transformation method together with (29) to initialize the Gauss–Newton method for pose estimation.
To evaluate the pose estimates, the position error is expressed by the Euclidean distance between the position vector r and its estimate r ^ , i.e., r e = r ^ r , while the orientation error is expressed by the axis-angle vector between the orientation matrix R and its estimate R ^ , i.e., u e = ( log ( R ^ R T ) ) [27]. Similar to the pose estimator derived from the calibration algorithm, we can derive the theoretical lower bounds for the pose estimate, i.e., MCRB ( T i ) and CRB ( T i ) , from (35) and (36), respectively, by restricting the parameters Θ to T i . Then the left upper and right bottom 3 × 3 submatrices of MCRB ( T i ) are bounds on position and orientation, respectively. The same result applies to CRB ( T i ) .

4. Numerical Assessment

In this section, we first verify the performance of the proposed calibration algorithm and compare the mean squared error with the theoretical bound based on simulations. Then, the calibration procedure is used to determine the intrinsic parameters in a simulation setup, after which the pose is estimated. To obtain a comparative study, this section gives a comparison between our proposed pose estimation method and the camera’s PnP and iterative methods.

4.1. Calibration

As the number of parameters to be estimated in the calibration phase is large, to have good performance, the number and quality of observations must be sufficiently high. To ensure that these concerns are met, we consider a dedicated calibration setup, where we place the LEDs in a plane at vertical distance sufficiently close to the receiver, i.e., z = 0.65 m, to obtain a sufficiently high SNR. Furthermore, to guarantee the receiver has sufficient LEDs within its FOV, we assume N L = N L , r 2 LEDs are placed in a square grid with an area of 400 cm 2 , where the LED in the i th row and j th column has position r L = 20 N L , r 1 i 1 j 1 0 10 10 0 cm, i , j { 1 , , N L , r } . All LEDs are assumed to have a transmitted power P t = 1 W, and a Lambertian pattern with order γ = 1 . We set the origin of the system frame as the center of the LED plane.
The receiver consists of a circular QPD with radius r p = 5 mm, resulting in an active area A p = π r p 2 , and responsivity R p = 0.4 A/W. Above the QPD, we place at a height | h | = 3.0 mm a circular aperture with radius l = 2.5 mm. In our simulations, we set the misalignment between the centers of the aperture and the QPD equal to k m , q = 0.5 , 0.3 mm. The noise in the different quadrants of the QPD is assumed to have a correlation coefficient ρ = 0.7 . The orientation R of the receiver is defined with the ZXZ Euler angles [ θ α , θ β , θ γ ] [31]. We assume the azimuth and roll angles, i.e., θ α and θ γ , are uniformly distributed over [ 0 , 2 π ) , and we set the elevation angle θ β as a variable. In this way, the receiver has a random orientation but a controlled tilt angle of θ β . After fixing the orientation, we still need to specify the position of the receiver. To ensure that most LEDs are observed, we select the position of the receiver so that it is pointing toward the center of the LEDs plane. This corresponds to the position r = z R e 3 e 3 T R e 3 .
In the simulations, we use (15) to generate the observations. We assume the shot noise has power spectral density N 0 = 2.10 × 10 22 A 2 / Hz , which corresponds to a background spectral irradiance p n = 5.8 × 10 6 W / ( cm 2 · nm ) and a visible light bandwidth Δ λ = 360 nm [32]. Supposing the electrical bandwidth equals B = 1 MHz [8], we can compute the noise variance σ w 3 2 with σ w 3 2 = N 0 B . Simulation results will be given in terms of the effective SNR = ( γ + 1 ) A a R t P t 2 π σ w 3 2 with A a = π l 2 .(As the received SNR depends on the receiver pose, a fixed received SNR will limit the parameter space of the receiver pose. As this will complicate the simulation, we chose to fix the effective SNR) Given the optical transmitted power P t = 1 W, we obtain the effective SNR = 44.73 dB. In our simulations, we will use the range SNR [ 25 , 65 ] dB to take into account variations of the system parameters.
First, we evaluate the performance of the designed calibration estimator as a function of SNR. As a baseline method, we consider the calibration algorithm for cameras [33]. Because the algorithm of [33] takes the image of an object as input, which in this case corresponds to the observation of the light spot’s position, we use (29) prior to this algorithm to convert the observation t to the estimated center of the light spot. The number of LEDs in the LED plane is set to N L = 25 , and the receiver observes the plane from N T = 4 randomly generated poses with θ β = π / 9 rad. We vary the SNR from 25 to 65 dB and plot in Figure 5a the resulting root mean square error (RMSE) (denoted by CAL) versus ξ ¯ , i.e., the received SNR averaged over the four poses. The root of the MCRB and CRB (denoted by rMCRB and rCRB) from Section 3.2 and the performance of the baseline method (denoted by BAS) are also plotted for comparison. As we can see from Figure 5a, at high ξ ¯ , the two calibration methods perform close to optimal, while at low ξ ¯ , the proposed calibration algorithm outperforms the baseline method. This can be explained as the bias in (29) for the baseline method is larger than for the proposed calibration method when ξ ¯ is low. The gap between MCRB and CRB reflects the net performance cost of the followed strategy, where we use the normalized differences instead of the RSS values and the approximated PDF. The gap shows that the robustness to imperfect knowledge of the transmitted power and radiation pattern of the LEDs, obtained by this strategy, is achieved at the cost of accuracy. The performance loss for the estimation of h is larger than for the estimation of k m .
Next, we evaluate the effect of the number N L of LEDs and the number N T of poses. We set SNR = 45 dB and θ β = π / 9 rad. Both variables influence the number of observations, so we expect that increasing either N L or N T will improve the performance of the calibration. In Figure 5b, we show the performance for N L , r = N L = { 2 , 3 , , 9 } with N T = 3 , and in Figure 5c for N T = { 2 , 4 , , 16 } with N L = 25 . The figures show that the lower bound is reached for both the proposed calibration method and the baseline method when N L or N T is sufficiently large. In that case, the performance improves logarithmically with N L and N T . When N L or N T are small, i.e., when the number of observations is small, the proposed method outperforms the baseline method.
Finally, we evaluate the effect of the orientation θ β of the receiver with respect to the LED plane. We set SNR = 45 dB, N L = 25 , and N T = 4 . The calibration performance for θ β varying from 0 rad to 5 π / 18 rad is shown in Figure 5d. When θ β is small, the performance first improves by increasing θ β . This can be explained as for small θ β , the randomly generated poses are very close to each other implying the solution set is tightly coupled and thus more susceptible to noise. For large θ β , the performance starts to degrade, as in our assumptions we assumed that the receiver is pointing to the center of the LEDs plane, implying large θ β corresponds to placing the receiver far from the center of the LEDs plane. As a consequence, the RSS values will experience lower SNR values, causing the slight degradation of the performance. The optimal performance will therefore occur at intermediate values of θ β .

4.2. Pose Estimation

In the pose estimation phase, we assume the intrinsic parameters are known, i.e., in practice, we use the intrinsic values obtained in the calibration phase. As in this case, we only need to estimate the six parameters of the pose; the number of LEDs can be lower than for the calibration phase. Furthermore, we assume the VLP system is combined with the illumination system, so the number and placement of the LEDs must satisfy the illumination standards. To this end, we consider in the evaluation of the pose estimator a 5 m × 5 m × 3 m area, where N L = 25 LEDs with optical power P t = 5 W are mounted at the ceiling of the area, as shown in Figure 6, unless specified otherwise (assuming the optical power of the LEDs equals 5 W, an average horizontal illuminance of 500 lm (corresponding to an office area) in a plane at z = 2 m below the LEDs requires roughly 25 LEDs). We define the boundary vector b = 5 , 5 , 3 T , the number of LED columns in the X direction N L , X = 5 and Y direction N L , Y = 5 , and the positions of the LEDs are given by b 1 2 i 1 2 N L , X , b 2 2 j 1 2 N L , Y , b 3 T , with i { 1 , , N L , X } and j { 1 , , N L , Y } , where [ b ] i denotes the i th component of b . We assume the algorithm only uses the LEDs that result in all quadrants to produce a sufficiently large RSS output to ensure that the information in the RSS values is reliable. To take into account the effect of calibration and the possible uncertainty in the estimated parameters, we assume that the receiver has been calibrated and that the estimated intrinsic parameter vector k ^ N ( k , Σ k ^ ) , with Σ k ^ = diag ( [ σ u v 2 , σ u v 2 , σ h 2 ] ) , and σ u v and σ h are set to the RMSE of the calibration given in the Section 4.1.
To evaluate the performance of the proposed estimators, we consider a path with a circular pattern in the X Y plane and a sinusoidal pattern in the Z direction as shown in Figure 6. The radius of the circle in the X Y plane is 1.5 m, and the amplitude of the sinus is 0.2 m. The circle is centered at 2.5 , 2.5 , 1.2 T m. Starting at the coordinates 4.0 , 2.5 , 1.2 T m, the path oscillates sinusoidally in the Z direction and completes the path in three periods.
In Figure 7, we evaluate the proposed estimator with respect to the transmitted power. In this simulation, we set Σ k ^ equal to the corresponding RMSE values from the calibration when SNR = 55 dB, i.e., σ u v = 1.17 × 10 5 m and σ h = 2.60 × 10 5 m. The RMSE of the pose estimate, averaged over the path, is first evaluated in terms of the transmitted power P t = ( 1 + 2 m ) × 10 n W, where the variables m { 0 , , 4 } and n { 0 , 1 } specify the selected power values in the figure. Figure 7 shows this averaged RMSE of the designed estimator based on the least squares method (LSM) (denoted by LSM) and the root of the MCRB (denoted by rMCRB) as a function of the mean received SNR ξ ¯ . In comparison, we also add the performance of the proposed estimator when perfect knowledge of the intrinsic parameters is given (denoted by LSM-K) and when the intrinsic parameters are not calibrated and assumed to be the design values | h | = 2.8 mm and k m , q = 0 mm (denoted by LSM-X), while the true values are randomly selected from a Gaussian distribution with as averages the design values and variances given above. The averaged RMSE of the baseline method [33] (denoted by BAS) and the root of the CRB (denoted by rCRB) are also plotted for comparison. As a second baseline method, we also consider the classical PnP method (denoted by PnP). As expected, the performance of the estimator with perfect knowledge of the intrinsic parameters reaches the theoretical lower bound MRCB and improves with the SNR. On the other hand, the performance of the estimators without perfect knowledge of the intrinsic parameters, i.e., calibrated or not, shows an error floor when the received SNR increases. This error floor is caused by the uncertainty in the estimated intrinsic parameters, which becomes dominant when the received SNR is sufficiently large. Obviously, the error floor is much higher in the absence of calibration. Our algorithm strongly outperforms classical PnP, as expected, while the its performance is slightly better than the baseline method BAS. The gap between MCRB and CRB implies that the robustness to imperfect knowledge of the transmitted power and the radiation pattern of the LEDs, obtained by the strategy of normalizing differences, is achieved at the cost of accuracy.
To illustrate the robustness of the proposed algorithm to fluctuations in the transmitted power, we consider the situation where the transmitted power of an LED follows a Gaussian distribution [34], i.e., P t = κ P t , m , with κ N ( 1 , σ κ 2 ) . Experiments in [34] show that the standard deviation σ κ in general is smaller than σ κ < 0.0667 . In Figure 8, we plot the RMSE averaged over the entire path as a function of the standard deviation σ κ [ 0 , 0.0667 ] , assuming P t , m = 5 W, σ u v = 1.17 × 10 5 m, and σ h = 2.60 × 10 5 m. We observe that for the given range of σ κ , the RMSE is (essentially) independent of σ κ , implying the approach of the normalized differences is robust to fluctuations in the transmitted power.
Next, we show the RMSE of the pose estimates of the proposed estimator as a function of the uncertainty in k ^ for P t = 5 W and compare the resulting RMSE with that of the baseline method and with the theoretical lower bounds. Figure 9 shows that, as expected, the performance of the estimator with the perfect knowledge of the intrinsic parameters reaches the theoretical lower bound and is independent of the uncertainty. The figure also illustrates that a large uncertainty degrades the performance of the estimators significantly. When the uncertainty becomes smaller, the performance of the estimators of LSM and BAS gradually reach and are bounded by the theoretical lower bound. This is expected, as due to the uncertainty, the estimators are in the best case as good as the ones with perfect knowledge of the intrinsic parameters. (Close to) optimal performance is achieved as long as the uncertainty is sufficiently small.
Finally, we evaluate the performance of the proposed receiver over the whole considered area. For the LEDs, we consider the setup of Figure 6, where the N L = 25 LEDs are attached in a square grid to the ceiling of the 5 m × 5 m area. We assume the receiver is placed at the positions belonging to a horizontal grid with a horizontal spacing of 0.2 m, at a vertical distance z d { 1 , 1.5 , 2 , 2.5 , 3 } m from the LEDs. Taking into account that the receiver in practice is pointing roughly upwards, we consider for each sample point a random orientation with an elevation angle θ β that is uniformly distributed within the interval 0 , θ β , m rad. We consider two scenarios, i.e., θ β , m = 0 rad for the case where the receiver is always pointing upwards, and θ β , m = π / 3 rad for the tilted scenario. Figure 10 shows the cumulative distribution (CDF) of the pose error when z d = 2 m. We observe that the proposed algorithm outperforms the baseline methods BAS and PnP. The figure further reveals that our algorithm converges approximately for 95% of the sample positions when θ β , m = 0 rad, i.e., when the receiver points upwards, while this reduces to 75% when a random tilted angle is considered, i.e., when θ β , m = π / 3 rad. Hence, as expected, the coverage degrades if the receiver has a random tilt. To evaluate the effect of the vertical distance z d on the performance, we compute the probability that the error is smaller than 1 m for the position and 0.2 rad for the orientation. Figure 11 shows that the performance of our estimator improves when the distance increases, as well as for the baseline method BAS. Despite the larger distance, and thus lower SNR, increasing the distance will bring more LEDs into the field-of-view of the receiver, implying more information is available to estimate the pose. On the other hand, for the PnP method, the performance degrades when the distance increases. This is explained as the PnP method is more sensitive to a reduction of the SNR, resulting in a performance degradation despite the larger number of LEDs available within the field-of-view.

5. Conclusions

In this paper, we investigate the pose estimation and calibration problem for the QADA receiver. To this end, the channel link is first modeled in terms of the receiver’s pose and misalignment using the perspective projection model. Then, we study the effect of the strategy of differencing and normalization to the probability density function (PDF), and it is shown that the applied strategy results in a complex PDF. Therefore, we derive an approximated PDF in a simple closed-form, based on which the calibration and the pose estimation algorithms using the least squares principle are proposed. The proposed calibration does not require any information about the pose of the receiver and is robust to the received power variation and imperfect knowledge of the LEDs’ radiation patterns, while the proposed pose estimator is able to directly estimate the pose from the (normalized) RSS values. Computer simulations confirm that the proposed calibration algorithm is effective and that the proposed pose estimator outperforms the PnP algorithm.

Author Contributions

S.S. led the paper, drafted the main parts, derived the estimation methods; J.M.M.S. has supported the work on verification of the estimation methods; S.L. and H.S. were supervisors of the research and acquired the funding; H.S. proofread and edited the paper. All authors have read and agreed to the published version of the manuscript.

Funding

This work was partially funded by the EOS grant 30452698 from the Belgian Research Councils FWO and FNRS, by the Flemish Government (AI Research Program), by the National Natural Science Foundation of China under Grant 61771474, by the National Secretariat of Higher Education, Science, Technology and Innovation of Ecuador (SENESCYT).

Institutional Review Board Statement

Not applicable.

Informed Consent Statement

Not applicable.

Data Availability Statement

The study did not report any data.

Conflicts of Interest

The authors declare no conflict of interest.

Appendix A. The Expressions for the Gradient ∇Θμi,j

Based on the chain rule, Θ μ i , j = r S , q , i , j μ i , j Θ r S , q , i , j R 2 × 6 N T + 3 , where r S , q , i , j μ i , j is the gradient of μ i , j with respect to r S , q , i , j , and Θ r S , q , i , j is the gradient of r S , q , i , j with respect to Θ . r S , q , i , j μ i , j is determined by its k t h component
r S , q , i , j μ i , j k = 2 2 cos 2 α k e k T l π 1 e k T r S , q , i , j / l 2 ,
while Θ r S , q , i , j = k r S , q , i , j , T 1 r S , q , i , j , , T N T r S , q , i , j R 2 × 6 N T + 3 is given by
k r S , q , i , j = T i r ¯ L , j © e 3 T T i r ¯ L , j
and
T k r S , q , i , j = K e 3 T T i r ¯ L , j K T i r ¯ L , j e 3 T e 3 T T i r ¯ L , j 2 T i r ¯ L , j if k = i 0 R 2 × 6 if k i ,
where the gradient of r S , q , i , j with respect to T k (A3) is calculated with the infinitesimal perturbation, and the operators · © and · are defined as
a 1 , a 2 , a 3 , a 4 T © = a 3 0 a 1 0 a 3 a 2 ,
and
ξ T , η T = η I 3 ξ × 0 T 0 T ,
respectively.

References

  1. Hassan, N.U.; Naeem, A.; Pasha, M.A.; Jadoon, T.; Yuen, C. Indoor Positioning Using Visible LED Lights: A Survey. ACM Comput. Surv. 2015, 48, 1–32. [Google Scholar] [CrossRef]
  2. Zhang, X.; Gao, Q.; Gong, C.; Xu, Z. User grouping and power allocation for NOMA visible light communication multi-cell networks. IEEE Commun. Lett. 2016, 21, 777–780. [Google Scholar] [CrossRef]
  3. Alsaadi, F.E.; Elmirghani, J.M. High-speed spot diffusing mobile optical wireless system employing beam angle and power adaptation and imaging receivers. J. Light. Technol. 2010, 28, 2191–2206. [Google Scholar] [CrossRef]
  4. Eroglu, Y.S.; Anjinappa, C.K.; Guvenc, I.; Pala, N. Slow Beam Steering and NOMA for Indoor Multi-User Visible Light Communications. IEEE Trans. Mob. Comput. 2021, 20, 1627–1641. [Google Scholar] [CrossRef]
  5. Eroglu, Y.S.; Guvenc, I.; Sahin, A.; Pala, N.; Yuksel, M. Diversity combining and piezoelectric beam steering for multi-element VLC networks. In Proceedings of the 3rd Workshop on Visible Light Communication Systems, New York, NY, USA, 3–7 October 2016; pp. 25–30. [Google Scholar]
  6. Rahaim, M.B.; Morrison, J.; Little, T.D.C. Beam Control for Indoor FSO and Dynamic Dual-Use VLC Lighting Systems. J. Commun. Inf. Netw. 2017, 2, 11–27. [Google Scholar] [CrossRef]
  7. Beysens, J.; Wang, Q.; Van Den Abeele, M.; Pollin, S. BlendVLC: A Cell-free VLC Network Architecture Empowered by Beamspot Blending. In Proceedings of the IEEE Conference on Computer Communications, INFOCOM 2021, Vancouver, BC, Canada, 10–13 May 2021; pp. 1–10. [Google Scholar]
  8. Cincotta, S.; He, C.; Neild, A.; Armstrong, J. High angular resolution visible light positioning using a quadrant photodiode angular diversity aperture receiver (QADA). Opt. Express 2018, 26, 9230–9242. [Google Scholar] [CrossRef]
  9. Cincotta, S.; He, C.; Neild, A.; Armstrong, J. Indoor Visible Light Positioning: Overcoming the Practical Limitations of the Quadrant Angular Diversity Aperture Receiver (QADA) by Using the Two-Stage QADA-Plus Receiver. Sensors 2019, 19, 956. [Google Scholar] [CrossRef] [Green Version]
  10. Aparicio-Esteve, E.; Hernández, A.; Ureña, J.; Villadangos, J.M. Visible Light Positioning System based on a Quadrant Photodiode and Encoding Techniques. IEEE Trans. Instrum. Meas. 2019, 69, 5589–5603. [Google Scholar] [CrossRef]
  11. Aparicio-Esteve, E.; Hernández, A.; Ureña, J.; Villadangos, J.M.; Ciudad, F. Estimation of the Polar Angle in a 3D Infrared Indoor Positioning System based on a QADA receiver. In Proceedings of the 2019 International Conference on Indoor Positioning and Indoor Navigation (IPIN), Pisa, Italy, 30 September–3 October 2019; pp. 1–8. [Google Scholar]
  12. Aparicio-Esteve, E.; Hernández, A.; Ureña, J. Design, Calibration, and Evaluation of a Long-Range 3-D Infrared Positioning System Based on Encoding Techniques. IEEE Trans. Instrum. Meas. 2021, 70, 1–13. [Google Scholar] [CrossRef]
  13. Mohammed, M.; He, C.; Cincotta, S.; Neild, A.; Armstrong, J. Communication aspects of visible light positioning (VLP) systems using a quadrature angular diversity aperture (QADA) receiver. Sensors 2020, 20, 1977. [Google Scholar] [CrossRef] [Green Version]
  14. Aparicio-Esteve, E.; Ureña, J.; Hernández, A.; Pizarro, D.; Moltó, D. Using Perspective-n-Point Algorithms for a Local Positioning System Based on LEDs and a QADA Receiver. Sensors 2021, 21, 6537. [Google Scholar] [CrossRef] [PubMed]
  15. Lepetit, V.; Fua, P. Monocular Model-Based 3D Tracking of Rigid Objects; Now Publishers: Norwell, MA, USA, 2005. [Google Scholar]
  16. Lin, B.; Tang, X.; Ghassemlooy, Z.; Li, Y.; Zhang, S.; Wu, Y.; Li, H. An indoor VLC positioning system based on OFDMA. In Proceedings of the Asia Communications and Photonics Conference 2016, Wuhan, China, 2–5 November 2016. [Google Scholar]
  17. Lin, B.; Tang, X.; Ghassemlooy, Z.; Lin, C.; Li, Y. Experimental demonstration of an indoor VLC positioning system based on OFDMA. IEEE Photonics J. 2017, 9, 1–9. [Google Scholar] [CrossRef]
  18. Steendam, H.; Wang, T.Q.; Armstrong, J. Theoretical Lower Bound for Indoor Visible Light Positioning Using Received Signal Strength Measurements and an Aperture-Based Receiver. J. Light. Technol. 2017, 35, 309–319. [Google Scholar] [CrossRef] [Green Version]
  19. Hartley, R.; Zisserman, A. Multiple View Geometry in Computer Vision; Cambridge University Press: Cambridge, UK, 2003. [Google Scholar]
  20. Kahn, J.M.; Barry, J.R. Wireless infrared communications. Proc. IEEE 1997, 85, 265–298. [Google Scholar] [CrossRef] [Green Version]
  21. Marsaglia, G. Ratios of normal variables. J. Stat. Softw. 2006, 16, 1–10. [Google Scholar] [CrossRef]
  22. Curtiss, J.H. On the Distribution of the Quotient of Two Chance Variables. Ann. Math. Stat. 1941, 12, 409–421. [Google Scholar] [CrossRef]
  23. Hayya, J.; Armstrong, D.; Gressis, N. A note on the ratio of two normally distributed variables. Manag. Sci. 1975, 21, 1338–1341. [Google Scholar] [CrossRef] [Green Version]
  24. Zhang, Z. A flexible new technique for camera calibration. IEEE Trans. Pattern Anal. Mach. Intell. 2000, 22, 1330–1334. [Google Scholar] [CrossRef] [Green Version]
  25. Nguyen, H.T.; Kreinovich, V. Applications of Continuous Mathematics to Computer Science; Springer: Berlin/Heidelberg, Germany, 1997; Volume 38. [Google Scholar]
  26. Absil, P.A.; Mahony, R.; Sepulchre, R. Optimization Algorithms on Matrix Manifolds; Princeton University Press: Princeton, NJ, USA, 2009. [Google Scholar]
  27. Shen, S.; Li, S.; Steendam, H. Simultaneous Position and Orientation Estimation for Visible Light Systems with Multiple LEDs and Multiple PDs. IEEE J. Sel. Areas Commun. 2020, 38, 1866–1879. [Google Scholar] [CrossRef]
  28. Fortunati, S.; Gini, F.; Greco, M.S.; Richmond, C.D. Performance Bounds for Parameter Estimation under Misspecified Models: Fundamental Findings and Applications. IEEE Signal Process. Mag. 2017, 34, 142–157. [Google Scholar] [CrossRef] [Green Version]
  29. Rasmussen, C.E. Gaussian processes in machine learning. In Summer School on Machine Learning; Springer: Berlin/Heidelberg, Germany, 2003; pp. 63–71. [Google Scholar]
  30. Kay, S.M. Fundamentals of Statistical Signal Processing: Estimation Theory; Prentice-Hall: Upper Saddle River, NJ, USA, 1993. [Google Scholar]
  31. Pissanetzky, S. Rigid Body Kinematics and C++ Code; SciControls.com: Houston, TX, USA, 2005. [Google Scholar]
  32. Steendam, H. A 3-D Positioning Algorithm for AOA-Based VLP With an Aperture-Based Receiver. IEEE J. Sel. Areas Commun. 2018, 36, 23–33. [Google Scholar] [CrossRef]
  33. Bouguet, J.Y. Camera Calibration Toolbox for Matlab. 2008. Available online: http://www.vision.caltech.edu/bouguetj/calib_doc (accessed on 23 June 2022).
  34. Plets, D.; Bastiaens, S.; Martens, L.; Joseph, W.; Stevens, N. On the impact of LED power uncertainty on the accuracy of 2D and 3D Visible Light Positioning. Optik 2019, 195, 163027. [Google Scholar] [CrossRef]
Figure 1. Structure of QADA.
Figure 1. Structure of QADA.
Sensors 22 05073 g001
Figure 2. Geometry of the light spot of a light beam, where x S = [ r S , q ] 1 and y S = [ r S , q ] 2 .
Figure 2. Geometry of the light spot of a light beam, where x S = [ r S , q ] 1 and y S = [ r S , q ] 2 .
Sensors 22 05073 g002
Figure 3. The channel link model. The position of the LED and centers of the aperture and the light spot comply with the perspective projection model, while the shape and size of the light spot is determined by the aperture. Note that in this perspective projection model, the light spot does not correspond to the image of the LED, but it only relates the position of the light spot with the position of the LED.
Figure 3. The channel link model. The position of the LED and centers of the aperture and the light spot comply with the perspective projection model, while the shape and size of the light spot is determined by the aperture. Note that in this perspective projection model, the light spot does not correspond to the image of the LED, but it only relates the position of the light spot with the position of the LED.
Sensors 22 05073 g003
Figure 4. System frame for calibration.
Figure 4. System frame for calibration.
Sensors 22 05073 g004
Figure 5. RMSE of the estimator and the theoretical bound: (a) N L = 25 , N T = 4 , θ β = π / 9 rad; (b) SNR = 45 dB, N T = 3 , θ β = π / 9 rad; (c) SNR = 45 dB, N L = 25 , θ β = π / 9 rad; (d) SNR = 45 dB, N L = 25 , N T = 4 .
Figure 5. RMSE of the estimator and the theoretical bound: (a) N L = 25 , N T = 4 , θ β = π / 9 rad; (b) SNR = 45 dB, N T = 3 , θ β = π / 9 rad; (c) SNR = 45 dB, N L = 25 , θ β = π / 9 rad; (d) SNR = 45 dB, N L = 25 , N T = 4 .
Sensors 22 05073 g005
Figure 6. Simulation setup for the QADA receiver. The three orthonormal vectors in three different colors (the red, green, and blue vector represent the x-axis, y-axis, and z-axis, respectively) at each sample on the path represent the frame of the receiver. The pink arrows represent the LEDs (only a fraction of them are shown) on the ceiling. θ r indicates the traveled angle along the dotted circle in the X Y plane.
Figure 6. Simulation setup for the QADA receiver. The three orthonormal vectors in three different colors (the red, green, and blue vector represent the x-axis, y-axis, and z-axis, respectively) at each sample on the path represent the frame of the receiver. The pink arrows represent the LEDs (only a fraction of them are shown) on the ceiling. θ r indicates the traveled angle along the dotted circle in the X Y plane.
Sensors 22 05073 g006
Figure 7. RMSE for position and orientation estimates, compared with the mean received SNR ξ ¯ .
Figure 7. RMSE for position and orientation estimates, compared with the mean received SNR ξ ¯ .
Sensors 22 05073 g007
Figure 8. RMSE for position and orientation estimates as a function of σ κ .
Figure 8. RMSE for position and orientation estimates as a function of σ κ .
Sensors 22 05073 g008
Figure 9. RMSE for position and orientation estimates as a function of σ u v and σ h .
Figure 9. RMSE for position and orientation estimates as a function of σ u v and σ h .
Sensors 22 05073 g009
Figure 10. CDF of RMSE for position and orientation estimates for z d = 2 m: (a) position error; (b) orientation error.
Figure 10. CDF of RMSE for position and orientation estimates for z d = 2 m: (a) position error; (b) orientation error.
Sensors 22 05073 g010
Figure 11. Probability versus the distance z d for: (a) position error r e < 1 m; (b) orientation error u e < 0.2 rad.
Figure 11. Probability versus the distance z d for: (a) position error r e < 1 m; (b) orientation error u e < 0.2 rad.
Sensors 22 05073 g011
Publisher’s Note: MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Share and Cite

MDPI and ACS Style

Shen, S.; Menéndez Sánchez, J.M.; Li, S.; Steendam, H. Pose Estimation for Visible Light Systems Using a Quadrature Angular Diversity Aperture Receiver. Sensors 2022, 22, 5073. https://doi.org/10.3390/s22145073

AMA Style

Shen S, Menéndez Sánchez JM, Li S, Steendam H. Pose Estimation for Visible Light Systems Using a Quadrature Angular Diversity Aperture Receiver. Sensors. 2022; 22(14):5073. https://doi.org/10.3390/s22145073

Chicago/Turabian Style

Shen, Shengqiang, Jose Miguel Menéndez Sánchez, Shiyin Li, and Heidi Steendam. 2022. "Pose Estimation for Visible Light Systems Using a Quadrature Angular Diversity Aperture Receiver" Sensors 22, no. 14: 5073. https://doi.org/10.3390/s22145073

APA Style

Shen, S., Menéndez Sánchez, J. M., Li, S., & Steendam, H. (2022). Pose Estimation for Visible Light Systems Using a Quadrature Angular Diversity Aperture Receiver. Sensors, 22(14), 5073. https://doi.org/10.3390/s22145073

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Metrics

Back to TopTop