Next Article in Journal
Study on Spatial and Temporal Characteristics of Surface Albedo at the Northern Edge of the Badain Jaran Desert Based on C + STNLFFM Model
Next Article in Special Issue
Comparison of Different Cure Monitoring Techniques
Previous Article in Journal
Enabling Blockchain Services for IoE with Zk-Rollups
Previous Article in Special Issue
Microsphere Coupled Off-Core Fiber Sensor for Ultrasound Sensing
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

Single LED, Single PD-Based Adaptive Bayesian Tracking Method

Department of Electrical and Computer Engineering, Sungkyunkwan University, Suwon 16419, Korea
*
Authors to whom correspondence should be addressed.
Sensors 2022, 22(17), 6488; https://doi.org/10.3390/s22176488
Submission received: 28 July 2022 / Revised: 20 August 2022 / Accepted: 24 August 2022 / Published: 29 August 2022

Abstract

:
Recently, with the growing interest in indoor location-based services, visible light positioning (VLP) systems have been extensively studied owing to their advantages of low cost, high energy efficiency, and no electromagnetic interference. However, due to structural limitations which lead to the requirement of multiple signal sources, it has been challenging to apply VLP in real-world scenarios. In this study, we propose a single LED, single PD-based tracking system that solves these problems by applying a new Bayesian method that can effectively reduce the computational burden of particle filters. The method of evaluating particle reliability developed in this work adjusts the number of particles on the fly. Using the absolute position of the single LED source, the long-term cumulative error of the inertial measurement unit can be continuously corrected. In this regard, the applicability of the VLP system can be enhanced in places where the multiple luminescent signals are hard to consistently detect. The proposed system was verified through experiments in a classroom and a corridor, and the results show an average error of less than 11 cm at travel distances of 80 to 100 m.

1. Introduction

In recent years, extensive research has been conducted on indoor positioning, where satellite signals are typically interrupted and degraded. Indoor positioning systems (IPSs) are configured through various types of signals, such as infrared and Wi-Fi, which have problems including short transmission distances, multipath effects, low accuracy, various types of signal interference, and the high cost of building such systems. To overcome these challenges, several studies have been conducted on visible light communication (VLC)-based positioning systems. By using light-emitting diodes (LEDs) as signal sources, VLC-based systems feature energy efficiency, low cost, long lifetime, and significant advantages in terms of construction costs because they use the LED infrastructure already built into the use environment [1,2,3,4]. In VLP systems, the most preferred method for low cost and high performance systems is the photodiode (PD)-based intensity-modulation/direct-detection (IM/DD) VLP system.
Because existing optical positioning systems are based on trilateration, at least three pieces of LED-to-PD signal information are essential [1,5]. This fundamental structure causes VLP systems to have the following limitations in real-life applications.
  • All LED signals must satisfy the line-of-sight (LOS) condition based on the element’s field of view (FOV) to receive multiple LED-PD signal information [1].
    -
    In Figure 1a, the area in which an optical signal can be stably received from each LED is fairly large (yellow circles ), while the area in which multiple signal-based positioning can be performed is quite limited (green trapezoids). Even if one or more signals can be stably received in the remaining area, positioning them is structurally impossible (red trapezoids).
    -
    LOS may not be guaranteed in various indoor and outdoor environments due to unexpected obstacles, user shadows, and extreme situations such as element failure (Figure 1b). If positioning is performed based on a single optical signal, the resulting area of positioning can be very wide.
  • When the radiation and incident angles exceed the FOV of each element, the error in the received signal gain increases exponentially [5]. Even if the FOV is adjusted to handle this problem, a large FOV results in poor accuracy owing to ambient or reflective light [6], whereas a small FOV results in the PD not receiving LED signals correctly [7].
  • There can be issues with LED management and reliability. Depending on the device characteristics, the junction temperature of an LED may change due to the driving current, ambient temperature, and self-heating. The main wavelength, signal efficiency, and characteristics of the LED may change according to changes in temperature [8,9,10]. LED signal power may degrade over time, creating error directionality and adversely affecting overall positioning accuracy [1,11].
  • Additionally, the use of multiple LED channels induces interference effects such as crosstalk between channels [12,13,14] and inter-cell interference (ICI) [15,16,17].
  • Multiple PD-based VLP methods [18,19,20,21,22,23] have been developed to produce multiple LED-to-PD signals in a single-LED environment. These studies include a method using a horizontal PD arrangement spaced at arbitrary intervals [18,19], a method using angular diversity detectors [20,21], and a method combining the two methods [22,23]. However, these methods require bulky receiving devices, and it is difficult to secure the viewing angle for all PDs.
  • There exist a few recent studies exploiting a single LED signal. Image sensor arrays can capture the geometrical features of a single LED signal for complete identification [24,25]. However, such devices typically have limited scan rates (30–60 Hz) for the entire sensor array. It is difficult to increase the modulation frequency and extend signal bandwidths. This further limits the noise immunity of VLP systems. Li et al. employed a special LED lamp for VLP, although this required additional small luminescent beacons [26].
  • The general positioning environment, where the receiver and transmitter are not facing each other in parallel and the height is not limited, has not been sufficiently studied [2,6,27]. Because the optical signal pattern includes the cosine of the incidence angle ( ψ ) and the radiation angle ( ϕ ), the calculation becomes very complicated when the receiver is not horizontal. The effect of the X- and Y-axis movements is greater than the effect of the Z-axis movement, and it is relatively difficult to accurately estimate the height of the receiver. In our previous study [2], we investigated the existing problems with conventional VLP studies and analyzed in detail the reasons that these problems were difficult to solve. However, in [28], a new statistical model for device orientation was established based on the observation results of several participants using a smartphone. The authors showed that an average polar angle of 30 with a standard deviation of 9 occurred during walking activities.
Figure 1. Structural problem of general multiple signal-based VLP systems: (a) problems caused by restrictions on viewing angles and (b) problems caused by obstacles, failures, etc.
Figure 1. Structural problem of general multiple signal-based VLP systems: (a) problems caused by restrictions on viewing angles and (b) problems caused by obstacles, failures, etc.
Sensors 22 06488 g001
In general, inertial measurement units (IMUs) are among the most common elements found in electronic devices, and are mainly used to estimate the posture of the device. The inertial navigation system (INS) is a dead-reckoning navigation system that acquires dynamic information through direct measurement from an IMU sensor. INS is the most basic method of sensor tracking systems due to its simple implementation. However, because INS estimates the current position from the previous position with current measurements, there is a fatal problem when the starting position is unknown or the system is unexpectedly restarted. In addition, inertial measurement devices are vulnerable to drift and environmental errors in which the reference output value of the sensor is affected in a static state [29]. Although great performance can be shown in the short-term, the long-term error accumulates as the driving time increases by taking the form of integration [30].
In this study, we propose a single-LED, single-PD-based sensor tracking algorithm using an optical sensor and an IMU sensor. To the best of our knowledge, a method for three-dimensional (3D) position estimation exploiting a single optical signal has not been presented in any other study to date. The proposed method is a light-based adaptive particle filter (Li-APF) developed using a novel Bayesian tracking technique. A particle filter (PF) is a trial and error-based prediction algorithm that is highly effective for estimating the state of nonlinear systems with non-Gaussian noise distributions based on probability density distributions. However, the fundamental problem for precise positioning is that the number of calculations increases as the number of particles increases. The proposed method effectively reduces the computational cost of particle filters by utilizing optical signal characteristics. By combining particle positioning and INS, we implemented a complete tracking system that compensates for the drawbacks of each sensor. The characteristics of the proposed method are summarized as follows:
  • By requiring only one pair of optical signals, most of the fatal problems of conventional multi-source VLPs are solved. Because the proposed method requires the stable reception of only one signal, the operational coverage is very wide (Figure 1a, red trapezoids) and there are few restrictions due to obstacles (Figure 1b). Furthermore, by combining the INS results, positioning can be maintained even in a region where optical signals are blocked.
  • The proposed method can effectively adjust the initial distribution of the particle filter using the optical signal intensity information. We propose an index called “particle reliability” which evaluates the overall particle behavior. Through this index, the number of particles is flexibly adjusted; thus, the calculation burden of the particle filter is reduced. Online tracking is guaranteed for a nonlinear system with non-Gaussian noise by effectively reducing the computational cost of the particle filter.
  • Based on the absolute position of the LED signal, the system itself can estimate its own position even if it does not have any previous state information. The proposed algorithm corrects the long-term cumulative error of the INS based on the continuity of optical positioning.
The remainder of this paper is organized as follows. In Section 2, we describe the optical model used in the particle calculation process. Section 3 describes the proposed tracking algorithm. The overall flow is introduced and explained in detail for each part. In Section 4, we verify the particle update process, which is the core of the proposed algorithm, through simulation. In Section 5, we describe an indoor tracking experiment with the proposed method using a field-programmable gate array (FPGA) and simple hardware, then discuss its numerical results. Finally, we present our conclusions in Section 6.

2. Optical Model Representation in Cartesian Coordinate System

The proposed algorithm uses the ideal optical signal strength value according to the position and tilt of each particle in the initialization and particle updating processes. Figure 2 shows an optical signal transmission system. Assuming that the position of the LED is L O x L , y L , z L and the position of the device is X O x X , y X , z X , the corrected position of the device assuming the LED as the origin is P X , Y , Z = x X x L , y X y L , z X z L . The horizontal and vertical tilt of the PD is represented by T β , θ , and α is the horizontal angle at which the transmitter points toward the receiver, while ϕ and ψ are the radiation and incident angles with respect to the transmitter and receiver, respectively. When the unit vector of the receiver tilt and the direction of the device and the LED facing each other are A , B , and N , respectively, the incident angle ψ is expressed by the inner product of the two vectors as follows:
cos ψ = A · B A B = cos θ cos ϕ sin θ sin ϕ cos ( β α )
In addition, as P X , Y , Z = d × N , N = sin ϕ cos α , sin ϕ sin α , cos ϕ , the horizontal and vertical angles of the spherical coordinate system are expressed as follows:
cos ϕ = Z X 2 + Y 2 + Z 2 , sin ϕ = X 2 + Y 2 X 2 + Y 2 + Z 2 , cos α = X X 2 + Y 2 , sin α = Y X 2 + Y 2
Here, by applying the above to the Lambertian radiation pattern [31], the ideal optical signal strength H received from the device is derived as follows:
H ( P , T ) = ( m + 1 ) A 2 π · cos m ϕ cos n ψ · d 2 = K · Z m cos θ Z sin θ cos β X sin θ sin β Y n X 2 + Y 2 + Z 2 2 m n / 2 ,
where A is the physical area of the light detector and d is the distance between the LED and PD. In addition, m = ln 2 ln cos ϕ 1 / 2 and n = ln 2 ln cos ψ 1 / 2 , where ϕ 1 / 2 and ψ 1 / 2 are the half angles of the transmitter and receiver, respectively. By converting the signal strength model from a spherical coordinate system to a Cartesian coordinate system, which is a major factor in the particle update process, the calculation of signal strength can be simplified.

3. Single Visible Light Signal-Based Tracking Algorithm

In this study, we propose a new Bayesian tracking method, Li-APF. The proposed method can effectively reduce the computational cost of the particle filter by requiring a very small number of particles and managing them efficiently. Using the results of particle estimation and the INS, we compensate for the drawbacks of each result prior to the final decision. The proposed method uses the dead reckoning principle with inertial sensor data. By additionally referring to the absolute location information of the LED and the distance of the received signal strength (RSS), it is possible to operate without any initial information. Additionally, the INS error is compensated for by continuous positioning using the optical signal information.
Figure 3 shows an overall flowchart of the proposed tracking algorithm. The algorithm is divided into five steps. In the first step (①), raw sensor data are measured using an optical sensor and an IMU sensor, then the data are converted for use in an algorithm. Next, if the initial positioning is performed without any particle information, particle initialization (②) is performed; otherwise, the particle update process (③) is performed. The particle update process is the core process of the proposed algorithm, and includes particle unit movement, weight updating through the cost function, particle positioning, and resampling. The proposed method reduces the direct computational burden by efficiently reducing the total number of particles and further simplifies the calculation equation by modifying the optical signal intensity model to facilitate positioning analysis. In the next step (④), the particle positioning result is evaluated and the final position is determined according to the evaluation results. Finally, the INS error-correcting process (⑤) is conditionally added. Each of these five processes is described in detail in the following subsections. The primary parameters are listed in Table 1.

3.1. Sensor Data Acquisition and Analysis Process

First, environmental data are collected through PD and IMU sensors and the values are processed for use in the positioning algorithms. In the case of IM/DD-based optical signals, signals from various channels are combined and received; in general, the signal strength H l t for each channel is classified by applying multiplexing and demultiplexing. The reference signal ( H t ) that determines the device position was selected for the largest received optical signal intensity ( H l t ). The IMU sensor data output the local data, that is, the data based on the sensor orientation, meaning that it should be converted into global coordinates. To this end, the tilt of the sensor ( θ , β ) is first estimated through the Attitude and Heading Reference System (AHRS) based on the sensor data, then the sensor data are converted into global coordinates by reverse correction of the posture. In addition, the movement speed V INS and unit time movement distance S INS of the current device are calculated through the accumulation of acceleration data.

3.2. Particle Initialization Process

After analyzing the sensor data, the next step is to check whether there are previously generated particles, and if not, to initialize the particles. Examples of particle initialization with both the general PF method and the proposed method are shown in Figure 4. The purple ovals in the figure indicate the positions of the LEDs. In the proposed system, the ideal optical signal strength can be calculated based on the distance from the LED, meaning that the initial random particle is set for the area that meets the signal strength condition (red dot), not for the entire map area (green dot). This condition effectively limits the initial particles. The initial position of particle P i 0 is set to satisfy the following condition:
P i 0 = r i · cos ( 2 π i / N P ) , r i · sin ( 2 π i / N P ) , h i H ( P i 0 , T 0 ) H 0 < ϵ H
where r i and h i are the distance between the LED and the particle in the x–y plane and z plane, respectively. Each particle has state information ( P i 0 , T 0 ) and weight values. In the initial stage, the weights of all the particles are initialized to the same W init value. Because there is no connection with the previous time, the particle update process cannot be performed, and the process continues to the next acquisition.

3.3. Particle Updating Process

If previously generated particles exist, the particle updating process is performed. In this process, various probability particles distributed in space are virtually moved using the unit time movement information from the IMU sensor data. The probability weight of each particle is updated using the cost function calculated by comparing the ideal value with the measured value of the previous and current optical signal strengths. Thereafter, the weighted average of the updated particles is determined as the current position, and particles with low weight are discarded and regenerated randomly.
When the previous and current states of the i-th particle are ( P i t 1 , T t 1 ) and ( P i t , T t ) , they can be derived as P i t = P i t 1 + S INS . After moving all particles, the cost function for each particle is calculated. The desired value determined through the cost function is the degree of the positional error of the positioning target, which is related to the distance. Because the distance d is a value proportional to the 1 / 2 power of the signal strength, the ratio is calculated by matching the measured signal strength and the ideal signal strength proposed in this study to the distance level. In addition, because the probability of the current particle and the probability of the previous particle must be considered together in order to derive the exact cost of the particle moving at a continuous measurement, it is expressed in one equation by multiplying the signal intensity ratio of the two values. Then, the square root calculation is performed because the whole expression is again expressed as a value of the square level of the distance, and its reciprocal makes a smaller error have a larger. Here, the cost function value C i of the particle is calculated as follows:
C i = 1 / H ( P i t 1 , T t 1 ) H t 1 1 · H ( P i t , T t ) H t 1
After calculating the cost function of each particle, the weight is updated. The weight is calculated as the sum of the previous weight and the current cost function at a constant ratio, r W : 1 r W (Equation (6)). Here, if the sampling frequency is f s and the stabilization time is T stable , the weight update rate r W is adjusted according to Equation (7). The stabilization time is arbitrarily designated by the user, and refers to the time required for the measurement of the accumulated data to affect the result by more than 90%:
W i t = W i t 1 · r W + C i · ( 1 r W ) = W init · ( r W ) t + + C i · ( 1 r W )
( r W ) T stable · f s < 0.1 r W < 0 . 1 1 / ( T stable · f s )
In this study, values of f s = 10 Hz, T stable = 0.5 s, and r W 0.631 were applied.
After updating the weight of all particles, the result of the particle positioning X P t is calculated by obtaining the weighted average value of all particles based on the updated weights. However, particles are scattered in various directions with respect to the LED. When all particles are considered, the average coordinate is inevitably designated as a location close to the LED. Therefore, the final position is calculated by considering only the upper 50% weights of the particles with the highest values. In addition, the absolute position compensates for the adjusted relative position based on the LED position. When the set of particles with the upper 50% weights is U , the estimated position X P t is derived as follows:
X P t = W i U t * × P i U t + L O ,
where W i U t * represents a normalized weight of upper 50% particles.
Next, the particle reliability R P t is calculated using the weighted average of the cost function based on the updated weights of all particles and the currently received signal strength, as follows:
R P t = i N P W i t * C i · H t ,
where W i t * represents the normalized weight of all the particles. By taking the weighted average, the index of reliability can be defined using the average cost function and the degree of positioning error accordingly, further considering that the signal strength becomes smaller as the distance from the LED increases. Because the reliability is based on the average cost function, it is considered to be an indicator of whether the overall behavior of particles matches the measured signal information.
In this study, the expected standard reliability R typ and minimum reliability R min are explicitly defined, allowing the current particle reliability R p t to be evaluated during the Li-APF procedure. These parameters are used to determine the particle numbers and the resulting position from INS and PF estimates. When the arbitrary representative reference coordinates of the previous and current position are P typ t 1 , P typ t and the actual device location is P O t 1 , P O t , the ideal values of the measured signal intensity are expressed as H t 1 = H ( P O t 1 , T t 1 ) and H t = H ( P O t , T t ) . If the distance calculated by any P k is D k = ( X k 2 + Y k 2 + Z k 2 ) 1 / 2 , the representative cost function can be expressed as (10). This equation is constructed based on (5) mentioned above.
C typ = 1 / H ( P typ t 1 , T t 1 ) H ( P O t 1 , T t 1 ) 1 · H ( P typ t , T t ) H ( P O t , T t ) 1 = 1 / D typ t 1 D O t 1 D O t 1 · D typ t D O t D O t
Here, if the difference between the representative coordinates and the actual coordinates, that is, the error distance, is Δ typ t 1 = D O t 1 D typ t 1 , Δ typ t = D O t D typ t , and Δ typ t 1 Δ typ t = Δ typ , D O t 1 D O t = D O , then C typ can be further simplified as below. As the sampling times of sensor tracking systems are short enough, this assumption is convincing.
C typ = D O t 1 · D O t Δ typ t 1 · Δ typ t D O Δ typ
By setting the reference conditions as ϕ = ϕ 1 / 2 and ψ = ψ 1 / 2 and applying this to Equation (9), the standard reliability can be derived as follows:
R typ = C typ · H P O t , T t K · cos m ϕ 1 / 2 cos n ψ 1 / 2 · 1 Δ typ = K 2 · 1 Δ typ .
This equation shows that R typ is inversely proportional to the reference error that the user can arbitrarily designate. Similar to the representative reliability, the minimum reliability is derived as follows:
R min K 2 · 1 Δ min .
In this study, we apply Δ typ = 0.05 m and Δ min = 0.5 m for R typ and R min , respectively.
Next, we calculate the weighted standard deviation of the particle angle σ P , which can determine the degree of distribution of all particles. Unlike particle reliability, this value is an indicator that evaluates whether the estimated position X P t as a weighted average is sufficiently reliable. To calculate the degree of spreading of particles in an elliptical shape based on the LED, the angle between the estimated position and each particle, ω i , is calculated around the LED. The weighted standard deviation of the spread angle is then calculated.
ω i = atan ( X P t ) atan ( P i t )
σ P = i N W i t * ω i ω i ¯ 2
Here, atan ( · ) is a function representing an arctangent value on the x–y coordinate of the parameter vector based on the x-axis. In this study, we applied the minimum angle deviation as σ min = 45 .
Prior to particle resampling, the amount of computation can be efficiently reduced by adjusting the total number of particles, depending on the situation. The particle positioning performance can be classified into three sections based on the particle reliability value, namely, excellent, good, and bad. We set the minimum particle number that can guarantee accurate positioning when the particle performance is excellent. Conversely, we check the maximum number that allows the particle distribution to converge correctly within a certain time while the performance is bad. This study was designed to adjust the number of particles N P according to the following criteria:
N P = 10 , w h e n R P t > 0.8 R typ 30 , w h e n 0.8 R typ R P t > 1.2 R min 50 , w h e n 1.2 R min R P t
Through this process, if the number of particles decreases, the particles with the lowest weight are removed by that number, and if the number increases, an additional particle with W i t = 0 is generated. The position of the new particles is determined through the resampling process.
Particle resampling is the final step in the particle updating process, in which particles with low weights are discarded and new particles are created at random candidate locations. In this study, particles with W i t < 0.5 W init are rearranged to satisfy the following optical signal strength condition:
P i t = r · cos ( 2 π · R N ) , r · sin ( 2 π · R N ) , h H ( P i t , T t ) H t < E H ,
where R N is any random number between zero and one. The weights of all resampled particles are initialized to W init value.

3.4. Decision through Weighted Average with INS

The aforementioned particle updating process is considered to be based on the optical signal intensity information. However, when the reliability of optical signal reception is degraded due to defects in LED and PD elements or temporary signal distortion, it directly affects the entire system. To address this issue, we use an integrated INS tracking system. However, because the cumulative error is significant in the INS system, in this study the INS positioning result is determined by applying the current unit time movement from the previous integrated result, as follows:
X INS t = X Li APF t 1 + S INS
By deriving the positioning results of the particle filter and INS in each process and weighting them according to the ratio R P t / R typ : 1 , the final estimated position X Li APF t is determined as follows:
X Li APF t = R P t · X P t + R typ · X INS t R P t + R typ .
Figure 5 shows a brief portion of the operation of the proposed tracking system. In the figure, the purple triangle indicates the location of the LED. This expression is commonly applied to the following figures. The particle movement of the Li-APF and the current estimated state (particle-weighted mean) are shown as yellow lines and circles, and its reliability is represented as the size of the circle. The estimated states of the INS and weighted average decision are shown as blue and black circles. In the third estimation stage in the figure, because the received LED signal is sufficiently stable, the reliability of Li-APF is high, and the weighted average of the final estimation is focused on the Li-APF result. Additionally, owing to the high reliability, the total number of particles in the previous step decreases. Conversely, in the fifth estimation, the LED signal is relatively unstable, meaning that the reliability of Li-APF decreases, the number of particles increases, and the final judgment is biased on the INS result.

3.5. INS Calibration Using Particle Information

After estimating the current device’s position through the aforementioned processes, an additional process of calibrating the accumulated error of the INS is performed using particle information based on the optical signal. Figure 6 shows a problem that may arise due to the accumulation of INS errors, that is, the cumulative speed error. In this figure, it can be observed that the estimated result (black circle) is continuously derived according to the weighted average of the particle position (yellow circle) and the INS position (blue circle). In Figure 6a, which is largely influenced by the cumulative error of the INS, it can be seen that depending on the cost function, particle positioning operates in a standard manner for optical signal strength conditions, while the continuity of previous and current particles is not smooth enough (red arrow) due to the long-term error of INS. If this aspect continues, the reliability and deviation of the particles deteriorates and the error of the INS increases, resulting in fatal consequences for the overall system. Likewise, as the device moves away from the LED and the optical signal intensity weakens, the optical reliability decreases and the influence of the INS on the final positioning increases, leading to the same result. Figure 6b shows the state in which the INS error is corrected and the connection between the previous and current particles is improved. Even when the optical signal strength is weakened, this can guarantee overall positioning performance because the INS error has been corrected.
In this study, the difference between the consecutive points of the particles is used to compensate for this unit time movement error. For this, the previous particle positioning results and reliability should be stored in X P t 1 and R P t 1 . Then, the INS variable can be adjusted using the previous and current particle filter values, as follows:
V INS = V INS + λ INS · X P t S INS X P t 1
λ INS = λ ( R P t 1 · R P t / R typ 2 )
Here, the adjusting value is directly updated to V INS before being accumulated in the next data collection and analysis process. In this study, the acceleration correction ratio λ is applied as 0.1 s−1.
However, as can be seen from the flowchart in Figure 3, INS error correction uses the difference between the previous and current particle positioning results; thus, this process cannot be used when the reliability or distribution of the particle is poor.

4. Particle Update Simulation

In this section, we verify the functionality of the proposed algorithm in an environment where initial positional information is not specified. In this simulation, we focus on particle update behavior, not the entire process of the proposed algorithm. Therefore, this simulation result does not represent the final result of the proposed algorithm; the entire algorithm is practically evaluated in Section 5. We observe the movement of each particle and the convergence performance to the solution, and check the R P t , σ P , and X Li APF t values. The simulation is configured as follows. A single LED is located at a position (0 m, 0 m, 0 m) directly facing the floor. The device moves around the LED within a radius of 2 m, changing the vertical tilt randomly from 30 to + 30 with a speed of approximately 1 m/s. The height of the device is designated as 1.5 m. We measure the sensor data at a rate of 10 Hz for 3 s. In this simulation, an optical channel model with a signal-to-noise ratio of 30 dB is assumed and the noise from the IMU sensor is not considered. The representative error and the minimum error are set to Δ typ = 0.05 m and Δ min = 0.5 m, and thus the representative and minimum reliability values are R typ = 27.9 × 10 6 , R min = 278.7 × 10 6 . The reference deviation is set to σ min = 45 .
Figure 7 shows the particle updating process for each sampling time. In this figure and the following figure, the purple triangle indicates the location of a single LED. The green circle represents the position of each particle, and the size of the circle represents the cost function value of each particle proportionally. The numbers at the bottom of each particle represent integer digits of the percentage used for normalizing the updated weights of the entire particle. The red X represents the actual position of the moving device during the simulation. The X marked in yellow or blue is the position of the device calculated using the current particle information ( X P t ). Under the conditions of weighted deviation, if the particle positioning is judged to be valid with a small σ P σ min , it is marked in blue, and if it is not valid ( σ P > σ min ), it is marked in yellow. In the actual positioning process, invalid yellow X information is excluded. The current time T, particle reliability R P t , and weighted deviation σ P values are indicated at the top of each sample. In the initial state, fifty particles are uniformly spread to ideal device candidate positions according to the first received optical signal strength, all of which have weight values of the same size. Because there is no cost function calculation, the size of the circle is adjusted to have the same value. In addition, it can be observed that the initial particle distribution has an elliptical trajectory due to the tilt of the sensor ( θ = 29.32 , β = 230.45 ). After 0 s, the particle updating process begins. The cost function is large in particles near the actual device location, although particles in other locations show large cost functions as well. However, this pattern disappears after several particle updates. In addition, the number of particles is adjusted to thirty after 0.1 s, as the particle reliability R P t is greater than 1.2 R min . At 0.3 s and 0.9 s, the number of particles is adjusted to ten according to the R P > 0.8 R typ condition. Meanwhile, the value of each weighted deviation increases as the particles are evenly spread, and the deviation decreases as the particles become concentrated. In Figure 7, it can be seen that the distribution of particles after 1 s is relatively stable and there is little positioning error.
Figure 8 shows the particle initialization and position-determining process in various random environments. The particle updates for one simulation are accumulated and shown in each subfigure. In addition, for particles with information that is continuously updated without being removed at the previous time, a green line shows this continuity by connecting their previous and current positions. As shown in this figure, the actual position of the device can be accurately estimated within 1 s to at least 2 s using the proposed particle filter. It can be observed that the stable positioning result after 2 s has almost no error as compared the actual device position.

5. Experiment

In this section, we evaluate the tracking accuracy of the proposed algorithm through experiments in indoor environments. Although several LEDs can be used for experiments, our system consists of a single LED-based method that uses only the single strongest signal for each sample. In the experiments described here, the proposed system was validated in indoor environments with the user holding the sensor device and walking along a designated track with specific postures. By considering the mobile user statistics reported in [28], three postures determining the degree of sensor tilt and its variation were applied to each experiment. The experimental environment is shown in Figure 9. The same experiment was conducted in a classroom (14.4 m × 8.4 m × 2.5 m) and a corridor (17 m × 2.1 m × 2.8 m). To cover a large measurement space, multiple LEDs were used in each experiment; however, positioning was performed based on the strongest single signal for each sample time. In the classroom, eight LEDs were placed between (−3.6 m, −1.2 m) and (3.6 m, 1.2 m) at width and length intervals of 1.2 m. Measurements were made by moving three times counterclockwise at an average speed of 0.4 m/s along the outskirts of a rectangular area from (−4.5 m, −2.1 m) to (4.5 m, 2.1 m) and 1.2 m high. Six LEDs were placed in the corridor between (−6.0 m, 0 m) and (6.0 m, 0 m) at intervals of 1.2 m. The measurements were made by moving three times clockwise at an average speed of 0.4 m/s along a rectangular area from (−7.5 m, −0.8 m) to (7.5 m, 0.8 m) and 1.5 m high.
In this study, we implemented transmission and receiving modules that could be controlled separately by an FPGA. The transmission part included an LED board with a white LED (SJ-3W-CW) and a main board transmitting a power source and an FPGA signal, allowing several LEDs to be driven. The receiving board included a silicon PIN photodiode (OSD5-E) and an amplifier circuit, receiving an optical signal from the FPGA through a 12-bit A/D converter operating at 12 MHz. The optical signals were modulated and demodulated by the orthogonality of the Hadamard matrix, and we applied the Manchester code to filter out the effect of the ambient light source. The inertial sensor data were measured using an external IMU sensor (EBIMU-9DOFV4). Data from both sensors were obtained at intervals of 10 Hz. In addition, the free movement of the device was implemented by a self-synchronization operation through an asynchronous signal sequence in a wireless environment. The power of the receiving board was supplied through a portable battery. The receiver module is shown in Figure 10.
Data collection was carried out while holding and moving the device at the height of the pedestrian’s chest, that is, 1.3 m above the ground, to maintain it evenly. In addition, walking scenarios according to the pedestrian’s posture were divided into three categories. Walking postures were classified as walking with the device held horizontally (Pose-1, vertical tilt of 5 to + 5 ), walking while looking directly at the device (Pose-2, vertical tilt of + 25 to + 35 ), and walking while shaking the device freely (Pose-3, vertical tilt of 35 to + 35 ). The horizontal tilt was kept as stationary as possible based on the walking direction in the first two postures, and was changed randomly in free movement.
Figure 11 shows the differences in the measurement data according to the three walking postures. The figure shows part of the data measured in the corridor. On the left side of the figure, the black solid line, dotted line, and dashed line indicate the acceleration data in each X-, Y-, and Z-axis converted to global coordinates, and the brown line represents the vertical tilt of the device. The right side of the figure shows the change in intensity of the optical signal received from each of the six LEDs over time. This figure shows that the noise of the inertial data is small when there is little change in posture, that is, the tilt of the sensor does not change significantly (Poses 1 and 2). When there is a large change in posture (Pose 3), the noise in the inertial data is large. However, when the posture changes slightly, the change in the received optical signal intensity is small. When the change in posture is large, the change in the signal intensity is large as well. In a posture with a specific tilt (Pose 2), the overall optical signal intensity appears to be sloped to one side according to the incident angle of the optical signal received by the PD.
Figure 12 shows the results of the proposed algorithm in the classroom environment. The real trace made by connecting the four vertices is indicated by a black line. The yellow line represents the estimated position calculated during the particle update process, and the green line represents the position filtered according to the angle-weighted deviation condition ( σ P σ min ) among the results. The purple triangle indicates the location of a single LED. Here, the result removed by the reliability condition ( R P < 0.8 R min ) is not included. The blue line represents the final estimated position of the Li-APF algorithm determined using the INS and particles. As in the measurement data analysis, the experimental results show that the results of particle positioning are formed close to the actual track without significant errors in Pose 3, where there are many changes in signal strength; Poses 1 and 2, where signal changes are relatively limited, show a slight increase in error in the particle results. On the other hand, particle results with a large error can be observed among the results for Pose 3 that do not meet the angle deviation condition; these are caused by a relatively large INS error and a measurement error occurring during optical signal reception.
The proposed method configures each particle using P i t , P i t 1 and the moving distance. The cost function and weight are then determined by the current and previous sensor measurements ( H t , H t 1 ). Likewise, the next position of the particle ( P i t + 1 ) estimates the cost function and weight using H t + 1 and H t , and can be resampled according to Equation (17). Here, if the change in measurements is insignificant ( H t 1 H t H t + 1 ), the scattered individual particles show flat cost functions and particle resampling does not occur. Conversely, when the measurement changes steeply, the number of particles satisfying continuity conditions is reduced, and there are more chances to eliminate inappropriate particles from the solution during the resampling process. In this regard, it can be expected that the cases with Pose 3 should show the fewest positioning errors in the experiments. On the other hand, the errors of the INS may be increased due to the rapid change in acceleration. In light of all of these considerations, we can examine the synthesized positioning results according to different poses.
Figure 13 shows the tracking results in the corridor environment. Because the width of the space was relatively narrow, the distance between the LED and the trace was generally closer than in the classroom, and the corner was relatively short and changed direction quickly. Similar to the results in the classroom, the results with Pose 3 are stable, and the errors are relatively large in Poses 1 and 2. In addition, the results show that the deviation and positioning error of the particles are worse at the corner compared to the classroom.
Table 2 summarizes the errors of the INS, particle positioning, and Li-APF in each experimental environment. In this table, INS refers to the result with general inertial navigation positioning that uses only the IMU sensor data. In the position tracking system, the path error is calculated by defining the shortest error distance in all directions from the real trace. The table reports that INS is an inappropriate system for long-term positioning due to cumulative errors, and that the proposed particle positioning has an average error of 21.63 cm. Finally, the proposed Li-APF algorithm is reported to have excellent overall performance, with an average path error of 10.89 cm. Figure 14 shows the cumulative distribution function (CDF) of the Li-APF path error for each experiment. As previously observed, the error with Pose 3 is smallest in each place, and the result with Pose 1 with a lower INS error was better than that of Pose 2. In addition, in all cases it can be seen that the top 80% of results have an error within 10–20 cm.

6. Conclusions

In this study, a single LED, single PD-based sensor tracking algorithm was proposed with a novel light-based adaptive Bayesian tracking method. The proposed method effectively solved the structural limitations of the general VLP system by requiring only one signal pair and considering the tilt of the receiver. As the Li-APF borrows the concept of particle filters to effectively reduce the computational burden, the proposed system is robust for nonlinear systems with non-Gaussian noise distribution while being less burdensome on online tracking systems. The algorithm proposed in the was reviewed in depth via simulation. In addition, an experiment was conducted according to three device handling scenarios, moving a total of 79.2 m and 99.6 m in a classroom and corridor, respectively; the average error was found to be 10.891 cm. With the proposed method it is expected that real-life applications of VLP systems can be realized, with great potential in the development of various application fields.

Author Contributions

Conceptualization, D.K.; methodology, D.K.; validation, D.K., J.K.P. and J.T.K.; resources, D.K.; data curation, D.K.; writing—original draft preparation, D.K.; writing—review and editing, D.K., J.K.P. and J.T.K.; visualization, D.K.; supervision, J.K.P. and J.T.K.; project administration, J.K.P. and J.T.K. All authors have read and agreed to the published version of the manuscript.

Funding

This work was supported by a National Research Foundation of Korea (NRF) grant funded by the Korean government (MSIT) (No. 2022R1A2C1010879 and 2020R1A6A3A13075997).

Institutional Review Board Statement

Not applicable.

Informed Consent Statement

Not applicable.

Data Availability Statement

Not applicable.

Conflicts of Interest

The authors declare no conflict of interest.

References

  1. Kim, D.; Park, J.K.; Kim, J.T. High-Efficient and Low-Cost Biased Multilevel Modulation Technique for IM/DD-Based VLP Systems. IEEE Access 2020, 8, 218954–218965. [Google Scholar] [CrossRef]
  2. Kim, D.; Park, J.K.; Kim, J.T. Three-dimensional VLC positioning system model and method considering receiver tilt. IEEE Access 2019, 7, 132205–132216. [Google Scholar] [CrossRef]
  3. Chuang, Y.C.; Li, Z.Q.; Hsu, C.W.; Liu, Y.; Chow, C.W. Visible light communication and positioning using positioning cells and machine learning algorithms. Opt. Express 2019, 27, 16377–16383. [Google Scholar] [CrossRef] [PubMed]
  4. Cai, Y.; Guan, W.; Wu, Y.; Xie, C.; Chen, Y.; Fang, L. Indoor high precision three-dimensional positioning system based on visible light communication using particle swarm optimization. IEEE Photonics J. 2017, 9, 1–20. [Google Scholar] [CrossRef]
  5. 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]
  6. Kim, H.S.; Kim, D.R.; Yang, S.H.; Son, Y.H.; Han, S.K. An indoor visible light communication positioning system using a RF carrier allocation technique. J. Light. Technol. 2012, 31, 134–144. [Google Scholar] [CrossRef]
  7. Komine, T.; Nakagawa, M. Fundamental analysis for visible-light communication system using LED lights. IEEE Trans. Consum. Electron. 2004, 50, 100–107. [Google Scholar] [CrossRef]
  8. Lee, D.H.; Lee, S.Y.; Shim, J.I.; Seong, T.Y.; Amano, H. Effects of Current, Temperature, and Chip Size on the Performance of AlGaInP-Based Red Micro-Light-Emitting Diodes with Different Contact Schemes. ECS J. Solid State Sci. Technol. 2021, 10, 095001. [Google Scholar] [CrossRef]
  9. Ishii, R.; Koyama, Y.; Funato, M.; Kawakami, Y. Microscopic origin of thermal droop in blue-emitting InGaN/GaN quantum wells studied by temperature-dependent microphotoluminescence spectroscopy. Opt. Express 2021, 29, 22847–22854. [Google Scholar] [CrossRef]
  10. Mesleh, R.; Elgala, H.; Little, T.D. On the performance degradation of optical wireless OFDM communication systems due to changes in the LED junction temperature. In Proceedings of the International Conference on Telecommunications (ICT), Casablanca, Morocco, 6–8 May 2013; pp. 1–5. [Google Scholar]
  11. Chhajed, S.; Xi, Y.; Li, Y.L.; Gessmann, T.; Schubert, E.F. Influence of junction temperature on chromaticity and color-rendering properties of trichromatic white-light sources based on light-emitting diodes. J. Appl. Phys. 2005, 97, 54506. [Google Scholar] [CrossRef]
  12. Cui, L.; Tang, Y.; Jia, H.; Luo, J.; Gnade, B. Analysis of the multichannel WDM-VLC communication system. J. Light. Technol. 2016, 34, 5627–5634. [Google Scholar] [CrossRef]
  13. Zhang, L.; Wang, H.; Zhao, X.; Lu, F.; Zhao, X.; Shao, X. Experimental demonstration of a two-path parallel scheme for m-QAM-OFDM transmission through a turbulent-air-water channel in optical wireless communications. Opt. Express 2019, 27, 6672–6688. [Google Scholar] [CrossRef]
  14. Azhar, A.H.; Tran, T.; O’brien, D. A gigabit/s indoor wireless transmission using MIMO-OFDM visible-light communications. IEEE Photonics Technol. Lett. 2012, 25, 171–174. [Google Scholar] [CrossRef]
  15. Jung, S.Y.; Kwon, D.H.; Yang, S.H.; Han, S.K. Reduction of inter-cell interference in asynchronous multi-cellular VLC by using OFDMA-based cell partitioning. In Proceedings of the 18th International Conference on Transparent Optical Networks (ICTON), Trento, Italy, 10–14 July 2016; pp. 1–4. [Google Scholar]
  16. Zhou, K.; Gong, C.; Gao, Q.; Xu, Z. Inter-cell interference coordination for multi-color visible light communication networks. In Proceedings of the IEEE Global Conference on Signal and Information Processing (GlobalSIP), Washington, DC, USA, 7–9 December 2016; pp. 6–10. [Google Scholar]
  17. Maheepala, M.; Kouzani, A.Z.; Joordens, M.A. Light-based indoor positioning systems: A review. IEEE Sens. J. 2020, 20, 3971–3995. [Google Scholar] [CrossRef]
  18. Yang, S.; Jung, E.; Han, S. Indoor location estimation based on LED visible light communication using multiple optical receivers. IEEE Commun. Lett. 2013, 17, 1834–1837. [Google Scholar] [CrossRef]
  19. Steendam, H.; Wang, T.Q.; Armstring, J. Theoretical lower bound for indoor visible light positioning using received signal strength measurements and an aperture-based receiver. J. Light. Technol. 2016, 35, 309–319. [Google Scholar] [CrossRef]
  20. Yang, S.H.; Kim, H.S.; Son, Y.H.; Han, S.K. Three-dimensional visible light indoor localization using AOA and RSS with multiple optical receivers. J. Light. Technol. 2014, 32, 2480–2485. [Google Scholar] [CrossRef]
  21. Mmbaga, P.F.; Thompson, J.; Haas, H. Performance analysis of indoor diffuse VLC MIMO channels using angular diversity detectors. J. Light. Technol. 2015, 34, 1254–1266. [Google Scholar] [CrossRef]
  22. Yu, X.; Wang, J.; Lu, H. Single LED-based indoor positioning system using multiple photodetectors. IEEE Photonics J. 2018, 10, 1–8. [Google Scholar] [CrossRef]
  23. Han, W.; Wang, J.; Lu, H.; Chen, D. Visible light indoor positioning via an iterative algorithm based on an M5 model tree. Appl. Opt. 2020, 59, 10194–10200. [Google Scholar] [CrossRef]
  24. Cheng, H.; Xiao, C.; Ji, Y.; Ni, J.; Wang, T. A Single LED Visible Light Positioning System Based on Geometric Features and CMOS Camera. IEEE Photonics Technol. Lett. 2020, 32, 1097–1100. [Google Scholar] [CrossRef]
  25. Zhang, R.; Zhong, W.D.; Kemao, Q.; Zhang, S. A Single LED Positioning System Based on Circle Projection. IEEE Photonics J. 2017, 9, 7905209. [Google Scholar] [CrossRef]
  26. Li, H.; Huang, H.; Xu, Y.; Wei, Z.; Yuan, S.; Lin, P.; Chen, Z. A Fast and High-Accuracy Real-Time Visible Light Positioning System Based on Single LED Lamp with a Beacon. IEEE Photonics J. 2020, 12, 1–12. [Google Scholar] [CrossRef]
  27. Karunatilaka, D.; Zafar, F.; Kalavally, V.; Parthiban, R. LED based indoor visible light communications: State of the art. IEEE Commun. Surv. Tutor. 2015, 17, 1649–1678. [Google Scholar] [CrossRef]
  28. Soltani, M.D.; Purwita, A.A.; Zeng, Z.; Haas, H.; Safari, M. Modeling the random orientation of mobile devices: Measurement, analysis and LiFi use case. IEEE Trans. Commun. 2018, 67, 2157–2172. [Google Scholar] [CrossRef]
  29. Farrell, J.; Barth, M. The Global Positioning System and Inertial Navigation; Mcgraw-Hill: New York, NY, USA, 1999. [Google Scholar]
  30. Savage, P.G. Strapdown inertial navigation integration algorithm design part 1: Attitude algorithms. J. Guid. Control. Dyn. 1998, 21, 19–28. [Google Scholar] [CrossRef]
  31. Kahn, J.M.; Barry, J.R. Wireless infrared communications. Proc. IEEE 1997, 85, 265–298. [Google Scholar] [CrossRef] [Green Version]
Figure 2. Optical system representation in 3D space.
Figure 2. Optical system representation in 3D space.
Sensors 22 06488 g002
Figure 3. Overall flowchart of Li-APF algorithm.
Figure 3. Overall flowchart of Li-APF algorithm.
Sensors 22 06488 g003
Figure 4. Example of particle initialization (green dot: general PF; red dot: Li-APF).
Figure 4. Example of particle initialization (green dot: general PF; red dot: Li-APF).
Sensors 22 06488 g004
Figure 5. Brief representation of the process of the proposed algorithm.
Figure 5. Brief representation of the process of the proposed algorithm.
Sensors 22 06488 g005
Figure 6. Example of algorithm operation according to INS error: (a) when INS accumulated error occurs and (b) when the accumulated INS error is corrected.
Figure 6. Example of algorithm operation according to INS error: (a) when INS accumulated error occurs and (b) when the accumulated INS error is corrected.
Sensors 22 06488 g006
Figure 7. Timeline representation of the particle update process.
Figure 7. Timeline representation of the particle update process.
Sensors 22 06488 g007
Figure 8. Cumulative representation of particle updates in random trials.
Figure 8. Cumulative representation of particle updates in random trials.
Sensors 22 06488 g008
Figure 9. Experimental environment: (a) classroom (eight LEDs) and (b) corridor (six LEDs).
Figure 9. Experimental environment: (a) classroom (eight LEDs) and (b) corridor (six LEDs).
Sensors 22 06488 g009
Figure 10. Receiver module (PD, IMU sensor).
Figure 10. Receiver module (PD, IMU sensor).
Sensors 22 06488 g010
Figure 11. Part of the measurement data according to different postures.
Figure 11. Part of the measurement data according to different postures.
Sensors 22 06488 g011
Figure 12. Experimental results in classroom environment: (a) Pose 1; (b) Pose 2; (c) Pose 3.
Figure 12. Experimental results in classroom environment: (a) Pose 1; (b) Pose 2; (c) Pose 3.
Sensors 22 06488 g012
Figure 13. Experimental results in corridor environment: (a) Pose 1; (b) Pose 2; (c) Pose 3.
Figure 13. Experimental results in corridor environment: (a) Pose 1; (b) Pose 2; (c) Pose 3.
Sensors 22 06488 g013
Figure 14. Tracking error CDF of Li-APF algorithm in each environment.
Figure 14. Tracking error CDF of Li-APF algorithm in each environment.
Sensors 22 06488 g014
Table 1. Parameters of the proposed algorithm.
Table 1. Parameters of the proposed algorithm.
SymbolParameter
H l t 1 , H l t The intensity of the received optical signal of the l-th channel previous and current.
H t 1 , H t The reference (strongest) received optical signal strength previous and current.
V INS , S INS The instantaneous speed and unit time movement of the sensor.
θ , β Vertical and horizontal tilt of the device.
H ( · ) An ideal signal strength calculated according to a parameter.
H 0 , P 0 , T 0 The initial reference optical signal strength, position and tilt of the i-th particle.
W init Initialization value of particle weight.
P i t 1 , P i t The position of the i-th particle previous and current.
T t 1 , T t The tilt of the particles previous and current.
C i The cost function value of the i-th particle.
W i t 1 , W i t Weights of the i-th particles previous and current.
W i t * Normalized current weight.
R P t 1 , R P t Reliability of particles previous and current.
σ P Weighted deviation of the particle angles.
P O t 1 , P O t The virtual position of the previous and current device.
P typ t 1 , P typ t , P min t 1 , P min t The virtual reference position of the previous and current device.
Δ typ t 1 , Δ typ t , Δ typ The virtual position errors and generalization values of the previous
Δ min t 1 , Δ min t , Δ min and current reference device.
C typ , C min The reference cost function value calculated through the virtual device position.
R typ , R min The reference reliability calculated through the reference cost function value.
σ min The reference minimum weighted deviation of the particle angles.
N P Total number of particles.
X P t 1 , X P t The results of the particle positioning previous and current.
X INS t 1 , X INS t The result of the INS positioning previous and current.
X Li APF t 1 , X Li APF t The result of the Li-APF positioning previous and current.
Table 2. The results of the tracking experiment in various environments.
Table 2. The results of the tracking experiment in various environments.
EnvironmentAverage Error [cm]
(Place/Tilt)INSParticleLi-APF
ClassroomPose-1 ( 5 5 )122.52625.45412.515
Pose-2 ( 25 35 )365.74627.57814.367
Pose-3 ( 35 35 )527.97018.1418.569
CorridorPose-1 ( 5 5 )153.77519.21010.222
Pose-2 ( 25 35 )203.25822.89011.014
Pose-3 ( 35 35 )607.65517.3328.658
Average330.15521.63410.891
Publisher’s Note: MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Share and Cite

MDPI and ACS Style

Kim, D.; Park, J.K.; Kim, J.T. Single LED, Single PD-Based Adaptive Bayesian Tracking Method. Sensors 2022, 22, 6488. https://doi.org/10.3390/s22176488

AMA Style

Kim D, Park JK, Kim JT. Single LED, Single PD-Based Adaptive Bayesian Tracking Method. Sensors. 2022; 22(17):6488. https://doi.org/10.3390/s22176488

Chicago/Turabian Style

Kim, Duckyong, Jong Kang Park, and Jong Tae Kim. 2022. "Single LED, Single PD-Based Adaptive Bayesian Tracking Method" Sensors 22, no. 17: 6488. https://doi.org/10.3390/s22176488

APA Style

Kim, D., Park, J. K., & Kim, J. T. (2022). Single LED, Single PD-Based Adaptive Bayesian Tracking Method. Sensors, 22(17), 6488. https://doi.org/10.3390/s22176488

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