An Improved Method for Localization of Wireless Capsule Endoscope Using Direct Position Determination

A novel methodology for localization of wireless-capsule-endoscope (WCE) by using the direct-position-determination (DPD) method is proposed. A WCE enables the diagnose of gastrointestinal tract disorders in a non-invasive visual manner. Conventional methods proposed for WCE localization are not optimal in terms of positioning accuracy since they include two disjoint stages: i) estimation of initial parameters such as direction, time, or time-difference of arrival (DOA, TOA, TDOA), and ii) localization performed by intersecting the loci provided by initial parameters. Moreover, most of these methods can localize only one signal transmitter. In contrast, the considered DPD is a single-step method which processes data from all sensor elements simultaneously and can resolve several (theoretically one less than sensors number) co-channel signal sources in a medical ward. More precisely, the proposed method can concurrently localize a WCE and some beacons attached to the patient’s body, which can supersede the use of micro-electromechanical-systems (MEMS). In addition, the considered DPD can locate some WCEs for multiple patients at the same time. To do this, DPD seeks for the peaks of a spatial profit function (SPF) which is based on the statistical cumulants of the signals observed at the array sensors placed around the medical ward. Here, we propose the classical DPD (CDPD) as a coarse but fast, and the generalized DPD (GDPD) as a fine and high resolution method for an improved localization of WCE. Also, the proposed technique can be easily extended to any other indoor applications, where the dimensions of the whole monitoring area and receiving array are about some carrier wavelengths. The superiority of CDPD over DOA-based WCE localization technique, and GDPD over CDPD is verified through comprehensive numerical analysis. We have also drawn the Cramer-Rao lower bound (CRLB) as a benchmark for performance analysis of GDPD. Moreover, we introduced a tighter lower bound on the localization mean square error (MSE) achieved with multiple group arrays (MGLB), corresponding to CDPD and DOA-based methods, to show the superiority of GDPD over these techniques.


I. INTRODUCTION A. BACKGROUND AND MOTIVATIONS
During the recent decades, there has been a growing interest on the localization of wireless-capsule-endoscope (WCE) by using radio-location techniques in the general context of wireless-medical-telemetry-services (WMTS) [1]- [6]. From the estimation theory point of view, the conventional WCE The associate editor coordinating the review of this manuscript and approving it for publication was Mohammad Zia Ur Rahman . localization techniques are sub-optimum since they have two separate and disjoint stages: first, the estimation of some initial parameters such as direction, time, time-difference of arrival (DOA, TOA, TDOA), or the received signal strength (RSS), and second, the intersection of outcome loci from those parameters to localize the target. The direct position determination (DPD) is proposed more recently as a single-step optimal technique for multiple co-channel signal localization by jointly using all captured data from sensor elements [7]- [11]. DPD is a statistical approach which uses VOLUME 9, 2021 This work is licensed under a Creative Commons Attribution 4.0 License. For more information, see https://creativecommons.org/licenses/by/4.0/ a spatial profit function (SPF) whose peaks indicate the location of the signal transmitters. More precisely, the SPF is the result of a mathematical solution to the localization problem, which is based on the statistical properties of signals sensed at the sensor array, usually their covariance or higher order cumulant matrices [12]. These solutions are commonly referred to as beamformers [13]. Nevertheless, the high accuracy of DPD method is achieved at the expense of some extra hardware requirements such as synchronous sampling of signals observed at all sensor elements [7]- [10].

B. LITERATURE SURVEY
The positioning accuracy of WCE using received-signalstrength (RSS) technique is investigated in [1], [4]. A DOAbased WCE localization technique by using an unscented Kalman filter is studied in [2], [14] in a scenario where the patient can move freely in the medical ward while a micro-electromechanical-system (MEMS) is used to recognize the body location and orientation. A performance and complexity analysis for a DOA/TOA/TDOA based WCE radio-location system with circular arrays is presented in [3]. More recently, a WCE localization based on the measurement of the magnetic field [5], and a two-way TOA oriented WCE radio-location technique utilizing finiteelement-method simulations [6] are proposed, where both techniques use a free shaped array on the patient's belly. However, the aforementioned bi-step radio-location methods are neither theoretically optimum nor can resolve co-channel signals. To the best of our knowledge, no DPD-based WCE or other biomedical sensor-based radio localization is proposed so far [15]- [18]. On the other hand, there are few works in classical DPD (CDPD) literature regarding indoor applications including WMTS or more specifically WCE localization. In their pioneering work [7], Weiss et al. considered a tactical field application, and showed that DPD coincides exactly with the maximum-likelihood (ML) estimation. In [8], [10], DPD is employed using different beamformers such as minimum-variance-distortionless-response (MVDR) and multiple-signal-classification (MUSIC), respectively. A performance analysis of DPD with both MVDR and MUSIC beamformers is addressed in [9], [11]. Some DPD applications for wideband systems like radar and orthogonal frequency-division multiplexing (OFDM) is provided in [19], [20]. More specifically, a performance analysis of DPD using MUSIC algorithm is presented in [21].

C. CONTRIBUTIONS
As shown in FIGURE 1, we define two types of array structures which will be used throughout this paper. More precisely, FIGURE 1-(a) shows the multiple-group-array (MGA) and  depicts the single-group-array (SGA), where the former is commonly-used in DOA-based radio location techniques. The main contributions of this paper for localization of a WCE and/or some beacons attached to the patient as an instance of WMTS, are summarized as follows. • We propose the CDPD, which is an optimum localization technique as an alternative to the commonlyused DOA-based techniques, both using the same MGA structure. CDPD is a medium resolution technique since it ignores statistical relations between received signals from any two different sensor groups.
• We propose the generalized DPD (GDPD) as a high resolution localization technique compared to CDPD, using an SGA antenna structure with the capability of simultaneous sampling on all sensor elements. Although GDPD provides us with a very higher accuracy than CDPD, it may suffer from some spurious peaks in SPF while locating emitters with co-channel signals.
• We show that our previously introduced technique called dynamic-sensor-array-response (DSAR) [22] can be used to effectively reduce the mentioned spurious peaks of SPF in GDPD.
• We introduce a lower bound on localization mean square error (MSE) achieved by MGA antenna structures (MGLB), which is tighter than the Cramer-Rao lower bound (CRLB) that is essentially defined for SGA structures (Appendix A).
• We introduce two visual analysis tools based on CRLB (CRgram) and its statistical gradient (CRgrad) at each location of the monitoring area, which can be used for visual performance evaluation.
• We propose different solutions to combat and decrease the fading effect due to multi-path propagation.
• We investigate the wave refraction phenomenon inside the human body and proposed two refraction canceling method based on DPD. First, a pre-cancellation technique that considers the refraction effect throughout the DPD algorithm, and second, a post-cancellation method that removes the estimation bias caused by the refraction (Appendix B) after running the DPD algorithm.

D. PAPER STRUCTURE
The rest of this paper is organized as follows. In Section II, we present the problem formulation for the proposed GDPD and CDPD methods by considering the commonly-adopted SGA and MGA structures. In Section III, we express the beamformers (SPF makers) that are then used in DPD formalism. In Section V, we compare the localization techniques based on the MGA (for DOA-based and/or CDPD) and SGA (for GDPD) array structures, where we also introduce some visual tools for better arrangement of array elements.
In Section VII, we investigate the effect of the human body and the refraction phenomenon due to wave traveling through two different environments. Then, we address the necessary modifications in the system model parameters to limit the refraction problem. In Section VI, we use the proposed localization methods to reduce the fading effect on the DPD algorithm. Subsequently, in Section VIII, we present our simulations and numerical results to validate the superiority of our proposed methods in terms of positioning accuracy and resolution. Finally, we conclude the paper in Section IX.

E. NOMENCLATURE
The following conventions are used throughout this paper. Dirac delta function Throughout this paper, the DOA or azimuth angles are considered in navigation, (i.e., clock oriented and not trigonometric) mode.

II. PROBLEM FORMULATION
We consider an MGA with M separated L-element arrays/groups with element spacing about λ. For an indoor application such as a medical ward, the aforementioned MGA can also be considered as an SGA with U = LM sensor elements as depicted in FIGURE 1. We also consider K emitters transmitting uncorrelated co-channel signals. All receiver sensors are assumed to be stationary during the observation time interval, which is much larger than the maximum possible propagation delay over receiver sensors array. Meanwhile, the emitted signals are also assumed stationary during observation time since this condition can be met by human walking speed through a medical ward. The three dimensional (3D) coordinate vector of k-th transmitter is denoted by p k = [x k , y k , z k ] T where k = 1, . . . , K . We aim at estimating the unknown positions of K transmitters, based on the observations received by the sensors array. Notice that for a linear time-invariant (LTI) channel with impulse response h(t), the received x(t) corresponding to a narrow-band (NB) source s(t) = ae j2πf c t with a being the amplitude and f c the center frequency, writes which means that for a NB signal, the convolution reduces to a simple multiplication with channel response at signal's center frequency, which is a complex scalar.
Considering a U -element SGA, we model the channel response between the u-th receiver sensor and the k-th signal emitter as where f c is the carrier frequency and indicates the stochastic part of the channel response where 1 represents the line-of-site (LOS) received signal component, and β uk is a zero-mean circularly-symmetric complex Gaussian random variable characterizing the fading effect due to multi-path propagation for instance, and finally denotes the deterministic part of channel response, with g uk and δ uk as amplitude and phase responses of the channel respectively, which can be expressed in terms of the channel path-length d uk , and the normalization factors g 0 and d 0 . Therefore, the statistical NB signal model can be expressed as where x and b are vectors of random variables (RV) characterizing the received signal and the generated noise at U sensor elements, respectively, P = [p 1 , . . . , p K ] T is the coordinates matrix of K transmitters, and s is a vector of RVs denoting the transmitted signals from K emitters. We have with A, α, and A being the channel response, the fading factors, and the sensor-array-response or pointing matrices, respectively. The k-th column of array response matrix A that corresponds to the k-th emitter's location, is hereafter referred to as the pointing vector a k . Considering a slow fading or a quasi-static channel for this NB model, we can assume α to be constant during the signal observation time interval for each location estimation. VOLUME 9, 2021

A. CDPD FORMULATION
Considering an MGA structure that is mainly used in DOA-based localization techniques, we propose the CDPD as an improved alternative method that exploits MGA structure, where any statistical relations between observed data in different arrays are ignored. Consequently, the matrix form of (4) for the m-th array will be where all scalar RVs in (4) are replaced by their N -sample row vectors and so, RV column vectors x, s, and b in (4) are respectively replaced by matrices X m , S m , and B m . This leads to a least squares estimation (LSE) problem for M distinct arrays as where P is the estimated locations of K emitters and P is the vector of K scanning points through the monitoring area.

B. GDPD FORMULATION
To achieve a method with more resolution than CDPD in an indoor localization application, we propose the GDPD technique that jointly takes into account the statistical dependency of all observations on receiver elements. This means that GDPD considers the SGA structure for all sensors that is based on the generalized form of the (6) as where Regarding the statistical dependcy of observations of whole sensor elements in an SGA, the LSE problem will be which is a generalized form of (7).

III. LOCALIZATION ALGORITHMS
Obviously, in the estimation problems such as (7) or (10) the A and S are unknown because we know nothing about the location and number of signal sources. We only can analytically explain the pointing matrix A as a function of K transmitters' locations. This means that at least a 3K -dimensional search is necessary to find the proper P, providing that K is known. To solve such a problem, we first consider (10) for a single transmitter (K = 1) and by ignoring the fading effect, we havep where a is the pointing vector of p, and s is the sample vector of emitted signal. We are looking for an SPF which depends only on the available data X and the scanning point p. To solve such a problem, we can estimate s as a minimizer to the Frobenius norm term in problem (11) which is given bŷ where a H a is the norm of pointing vector of p that should be normalized during the search of monitoring area ( a = 1). From (11) and (12), we havê and after some manipulations we get By defining the sample covariance matrix of observations as which is the well-known conventional or Bartlett beamformer whose output value indicates the received signal power by the array. This was first introduced as a proper solution for DOA estimation of a single transmitter [13], which we use in DPD.
To localize multiple signal emitters, we can exploit many improved beamformers that are originally introduced for DOA estimation of multiple sources. Some of these beamformers are mentioned in next Section. It is worth mentioning that these beamformers usually need the covariance matrix of observations, i.e., R. Thus, we consider this matrix for both CDPD and GDPD cases.
Regarding (15), the general form of observations covariance matrix for the GDPD will be N BB H denote the covariance matrices of source signals and whole elements noise, respectively.
Using analogy, the covariance matrix of observations for m-th array in CDPD will be with R S , and R Bm = 1 N B m B m H denoting the covariance matrices of source signals and m-th array's noise, respectively. The general form of observations covariance matrix R for CDPD problem is defined as which ignores the statistical dependency of M separated arrays.

A. BEAMFORMERS
All beamformers (SPF-oriented algorithms) proposed the in context of co-channel DOA estimation problems can be exploited in DPD due to the similarity of their problem formulations, for instance, Bartlett [13], Capon or MVDR [23], [24], adaptive angular response (AAR) [25], thermal noise algorithm (TNA) [26], and MUSIC beamformer [27]. In contrast, non-SPF-oriented algorithms such as ROOT-MUSIC and ESPRIT which provide the estimated DOA in a direct mathematical form, cannot be used in DPD. Here we present the formulations of some commonly-used beamformers as where a = a(p) is the pointing vector of point p that scans throughout the monitoring area to find the peaks of the SPF, one-by-one. Considering the uniform definitions of R for CDPD and GDPD, the solutions are uniquely formulated aŝ where k = 1, . . . , K , andp k is the estimated position of k-th signal emitter.

IV. CONSIDERATIONS ABOUT CONVENTIONAL WCE LOCALIZATION SCENARIOS
In a WCE lab such as that depicted in FIGURE 2, the first step of conventional localization technique is determining the DOA of the WCE signal by using each column (group) of sensors array around the WMTS ward, independently [2]. The second step is the estimation of intersection points of loci presented by each column array to estimate the WCE position. This is aided by using MEMS beacons positions to find the WCE location in the patient's body frame. An unscented Kalman filter is then used to improve target tracking while the patient walks through the medical ward. In contrast, our proposed CDPD and GDPD techniques, whose flowchart is depicted in FIGURE 3, localize the WCE in a single-step search. It can also obviate the need of using MEMS components, which yields a unified localization hardware.
Regardless of the localization technique, there are some considerations about this WCE lab specifications [2], that should be taken into account.
1) The distance between patient and each column arrays is not so long compared to the length of that array (2.5 m) to let us consider the array response as a function of DOA only. Therefore, the DOA-based localization contains an inevitable error in this scenario.
2) The sensor elements spacing b, is a more important parameter than the length of array. The optimum value  depends on λ and on the array geometry [13], which means that the array length may not be constant in contrast to the scenario in [2] which uses 3, 5, 7, or 10 sensors per column array with a length of 2.5m.
3) The far-field condition is violated through the WCE localization scenario mentioned in [2] since λ 21 cm for the standard 1430 MHz carrier frequency, while the distance from the walking patient can be less than 10λ = 2.1m.

V. COMPARING LOCALIZATION METHODS BASED ON MGA AND SGA STRUCTURES
It is obvious that in a localization problem, we expect a higher resolution and accuracy while using an SGA structure (for GDPD) compared to using an MGA (for CDPD), since the mutual information between different groups (arrays) is ignored in the latter. On the other hand, although both DOA-based localization techniques and CDPD exploit MGA structures, the CDPD performs very closer to MGLB since it jointly uses the information from all groups together in a single stage [7]. Therefore, we do not deal with the details VOLUME 9, 2021  of DOA-based localization methods here, and only compare CDPD as the best MGA-oriented technique to GDPD. We also consider that the GDPD may suffer from unwanted spurious peaks on its SPF, and show that our previously proposed DSAR technique can efficiently reduce these spurs from the SPF of GDPD. FIGURE 4 shows the details of a transmitter localization, through which we can qualitatively compare the normalized SPFs of (a) CDPD that is smooth and low-resolution, (b) GDPD presenting higher resolution but with possible spurious peaks around the target, and (c) GDPD+DSAR showing an ultra-high-resolution surface and also a high-contrast between the true peak and spurs. The DSAR method considers the free-space-loss (FSL) through the elements of pointing matrix A to provide a more realistic channel model. Also, FIGURE 5 depicts the general 2D heat map of CDPD SPF whose peak determines the location of signal emitter.

B. AN ANALYTIC COMPARISON: THEORETICAL LOWER BOUNDS
To compare the proposed methods in terms of positioning accuracy, we derived the CRLB and MGLB for GDPD and CDPD problems respectively in Appendix A, and showed that MGLB is a tighter lower bound (LB) than CRLB. We know that the inverse of the Fisher Information Matrix (J = F −1 ) contains the lower limits for covariance between each pair of parameters ψ = [x, y, z] T as To compare the performance of CDPD and GDPD, we consider the MGLB and CRLB respectively, for up to three transmitters in some arbitrary points on a horizontal plane in WTMS ward, namely (0,0,0)m, (2,0.4,0)m, and (2,-0.4,0)m. Comparing these two plots, it can be easily figured out that: • The theoretical resolution of GDPD is about 15 dB better than that of CDPD, which is due to the fact that GDPD exploits the whole covariance matrix R elements, while CDPD ignores the mutual information between different groups.
• GDPD is robust to the presence of multiple signal transmitters and also to the location of evaluation point through the monitoring area, while CDPD does not show such a robustness. Notice that the theoretical limit for number of signal sources is one less than the number of sensor elements in a group [13], which is L for CDPD and U = LM for GDPD. The aforementioned characteristics encourage us to propose the combination of CDPD and GDPD for localization of multiple biomedical sensors in a WMTS environment. It is worth mentioning that although GDPD has a very higher resolution than CDPD, it requires a huge computational burden without prior usage of CDPD as a coarse position estimator.

C. VISUAL TOOLS FOR ARRAY ARRANGEMENT ANALYSIS: CRgram AND CRgrad
Here, we introduce some visual analysis tools for finding the proper array geometry, using the elements of matrix J. For  any arbitrary location in the monitoring area, the CRLB is determined by computing tr(J) which involves only main diagonal elements [J] ii = σ 2 ψ i . In fact, the set of CRLB values for all grid points in a 2D monitoring area can make a surface or an image to show the quality of estimation at that location. We refer to the heat-map of this surface as the array Cramer-Rao-gram (CRgram) (see FIGURE 8). According to the gray color-map, the darker color indicates a lower error. Although the CRLB can be calculated for a 3D area, such an image can only be achieved for a 2D cut of the monitoring area usually perpendicular to one of the Cartesian axes {x, y, z}.
To provide a 3D analysis tool, we propose a line-based plot whose lines show the magnitude and direction of gradient of localization error. This plot uses the non-diagonal elements of J in addition to the diagonal ones. We define u as the CRLB gradient unit vector as and g as the CRLB gradient vector as  whose magnitude shows the mutual dependency of error distribution on the Cartesian axes, where ρ is the Pearson correlation coefficient, defined as We denote the resulted view or picture as array CRgrad (see FIGURE 9). Moreover, the combination of CRgram and CRgrad can be exploited as a visual tool to assess the error associated to an array, such as that depicted in FIGURE 10.

VI. REDUCING THE EFFECT OF CHANNEL FADING
To reduce deterioration due to the multipath fading channel, the first step consists in performing some physical regularities in the medical ward e.g., by covering the medical equipment with radiation absorbing materials [2]. In this way, the LOS signal becomes stronger than the reflections leading to less localization errors. The second step can be accomplished using mathematical techniques to compensate the fading impact. In Section II, we modeled the fading effect as a random factor in the channel response matrix, inspired by the model proposed by Weiss et al. in their pioneering works [7]- [9]. We have where u = 1, . . . , U and k = 1, . . . , K , and α uk is a circular complex Gaussian random variable with mean one and standard deviation σ f (fading depth), drawn independently per u-th sensor and k-th emitter (i.e., per each channel path). considering a Rician fading channel for each path we have It is worth mentioning that we assume slow fading or quasi-static channels in such a way that each α uk can be treated as a constant during each observation time interval for one DPD procedure. We modeled the fading effect in our system model in Section II, by Hadamard multiplying the fading factor matrix α = [α uk ] to array response (pointing matrix) A. According to (29), it is clear that averaging can asymptotically neutralize the fading coefficients to a constant value. Assuming that DPD is an unbiased estimator, we propose three possible averaging-oriented schemes to overcome the fading effect: first, the averaging over coordinates (Avg_p), second, the averaging over pointing vector (Avg_a) of each estimated location (both with heavy computational loads), and third and most important one, the averaging over covariance matrix R (Avg_R) which we discuss as follows.
A. AVERAGING OVER COVARIANCE MATRIX R Considering (5) and (17), and assuming uncorrelated noise at receivers the sample covariance matrix of observations is where σ 2 is the noise power at receiver sensors. Providing that the source signals are uncorrelated we have R S = diag[P 1 , . . . , P K ] and the elements of R will be where u, v = 1, . . . , U and k = 1, . . . , K . Considering (29), the statistical expectation of (31) with respect to α uk will be where with σ f as fading depth, and therefore we have in which the fading depth σ f appears only on the last term. Therefore, we can rewrite (34) in a more compact form as where R is the fading-free covariance matrix we are looking for that equals the sum of first and second terms of (34), and P = diag([P 1 , . . . , P U ]), where P u = K k=1 |a uk | 2 P k equals the fading-free total received power at u-th sensor. Assuming the same P u = P for all receiver sensors, we have P = PI. Approximating the expected covariance matrix E{R} by averaging over R for N e consecutive estimation intervals, during each of which the α n matrix can be assumed constant (n = 1, . . . , N e ), we get We notify that the Eigen vectors of E{R} are the same as those of R but their Eigen values differ by σ 2 f P, which means that fading acts the same as a noise with power σ 2 f P. Moreover, notice that averaging over covariance matrix R is the fastest solution to reduce the fading effect compared to the two aforementioned techniques, since it needs only one pass of the DPD algorithm after computingR.

VII. THE EFFECTS OF THE HUMAN BODY ON WAVE PROPAGATION
Radio wave propagation through the human body leads to the refraction and also attenuation that change both phase and amplitude of the received signal at receiver sensors. Here, we focus on the refraction phenomenon.

A. REFRACTION
According to the Snell's law, changing the propagation environment causes the a traveling wave to be refracted. More precisely, the Snell's Law states that where n is the environment refraction index, θ indicates the incidence or refraction angle with respect to the normal vector In fact, the refraction phenomenon produces an intrinsic bias in a radio location problem, which means that if we don't regard the refraction in our localization algorithm, the signal source inside the human body can never be seen on its true place, even for a noise-free estimation. FIGURE 11 depicts a 2D schema of the virtual emitter displacement from coordinates p to p u , seen by the u-th sensor element, due the refraction phenomenon. Considering our discussion on a cylindrical model in Appendix B, and assuming an isotropic environment for human body, the final location resulting from averaging over displacement vectors corresponding to all U sensors is Therefore, the refraction leads to an overall intrinsic bias of Considering the array symmetry and a noise-free localization for a transmitter located at distance d from the cylinder axis, we have whered is the noise-free estimate of location. The bias value is a function of d as explained in Appendix B.

B. METHOD1: REFRACTION PRE-CANCELLATION
Here, we present a perfect solution for DPD of WCE while taking into account the refraction phenomenon that virtually changes the path-lengths for traveling waves. This solution is based on the modification of pointing vector (a) in (23) to compensate the refraction effect. FIGURE 12 depicts the wave refraction path using a cylindrical model for the human body. The signal is emitted from a hypothetical transmitter located at p, hits the cylinder surface at r u on angle θ i , leaves the surface on angle θ i , and finally reaches the sensor element at e u . The a r is the outwardpointing, unit radial vector of cylinder at r u . We indicate the plane containing {p, r u , e u } as the refraction plane.
To analyze the effect of the human body on localization problem, we remember that refraction affects both phase and amplitude responses of the channel between the transmitter (WCE) and the receiver sensor. In other words, refraction changes the elements of the pointing matrix A in (5) that plays an essential role in the localization problem. Notice that the wave propagation speed depends reciprocally on the refraction index of the environment, despite of the carrier frequency which is constant. Thus, λ which has a key role in characterizing the channel phase response, depends on the refraction index. In fact, this is the basis of the Snell's Law is where v stands for the wave propagation speed, and the subscripts have the same definition as those in (37). It is worth mentioning that the refraction plane contains the cylinder radius on point r u as a line segment. Considering the refraction path in FIGURE 12, we can determine the channel phase response between WCE and u-th sensor element aš where d ui = r u − p , d uo = e u − r u , andď u is the virtual distance between signal emitter at p and receiver element at e u . By comparingδ u to the free space phase response, we have To take the refraction effect into account, the elements of the pointing matrix A should be revised toǍ, by replacing δ u withδ u in (3). The coordinates of r u are determined by minimizing the cost function obtained from (37) by scanning r u on the cylinder surface. The cost function (44) is a more general form of the one presented in [2]. Unfortunately, there is no guaranty for the cost function (44) to have a valid minimum on the cylinder surface since the refraction prevents the emitted wave from any hypothetical point inside the cylinder to reach any arbitrary point outside. This means that some regions inside the cylinder are invisible to a sensor element. Fortunately, these regions have deterministic boundaries, which we investigate in Subsection VII-C.

C. HIDDEN AREAS: A CONSEQUENCE OF REFRACTION
Consider FIGURE 13, which shows a horizontal cut of the cylindrical model of human body. If a transmitter is located at the hashed regions, the sensor element at e u receives no signals from that transmitter, since the most lateral beam that reaches e u , is the one that hits the surface from the circle inside at critical angle (θ i = θ cr ) and leaves that tangentially (θ o = 90 • ). By definition, we have θ cr sin −1 ( 1 n i ), which is about 46 • for human body. In a more complex case, for a non-horizontal refraction plane with tilt angle η, we have where θ i is the horizontal projection of θ i . This means that other sensor elements in the same column as e u will have hidden areas with larger horizontal projection segments than that depicted in FIGURE 13. Obviously, this leads to deterioration of WCE localization performance due to the wrong DOA data grabbed from blind sensors.
To have a successful localization, it is necessary to determine all the sensor elements that cannot see the search point p (for all GDPD, CDPD, and DOA-based techniques), and then ignore those blind sensors observations to remove the irrelative data throughout the algorithm. To do this, a prior knowledge about the position and orientation of patient's body is necessary, which is gradually achieved during iterations of the algorithm mentioned in FIGURE 3. . An instance of hidden areas due to refraction on critical angle θ cr as an upper limit for inside incidence angle, which determines the hidden areas boundaries.

D. METHOD2: REFRACTION POST-CANCELLATION
In the tradeoff between the accuracy and simplicity, sometimes we prefer to have a lower computational load method even with lower performance. In fact, the pre-cancellation algorithm presented in Subsection VII-A consumes heavy computational resources. Therefore, we concentrate on the refraction post-cancellation idea, based on the removing of deterministic bias in (39). We analytically calculate the bias values with an acceptable resolution of d for only one time, and store them in a look-up-table (LUT). FIGURE 14 depicts the normalized value of the bias for an applicable range. Hence, to achieve a low computational load method, we localize the WCE ignoring the refraction effect, and then cancel the bias using the LUT. Notice that the bias LUT values are applicable as long asd is an invertible function of d. The aforementioned range (0 . . . 0.71R) seems realistic for the dimensions of intestine inside the body frame.

VIII. SIMULATIONS AND NUMERICAL RESULTS
Here, we investigate the improvement brought by using GDPD compared to CDPD and DOA-base localization techniques in a medical lab. The results of DOA-based methods can also be found in [1]- [3], [6]. We use the MUSIC beamformer for our simulations. The main parameter in our analysis is 3D location root-mean-square-error (RMSE) normalized to λ. Our results are averaged by using 100 Monte-Carlo experiments. There are 64 signal observations for each sensor element. The localization error is defined as the Euclidean distance between estimated locationp = (x,ŷ,ẑ) and known positions of target at p = (x, y, z) as We assume that the patient can walk through a determined area that guarantees the far-filed radiation conditions for both the WCE and beacons transmitting signals at carrier frequency 1430 MHz leading to wavelength as λ 21 cm. The standard frequency ranges 1395-1400 MHz or 1427-1432 MHz are approved by FCC Wireless Medical Telemetry Service (WMTS) [2], [28].
To show the vantage of exploiting our previously proposed DSAR technique [22], we consider both cases of having one and/or two transmitters (K = 1, 2). Notice that DPD methods are capable to locate co-channel signal sources. FIGURE 15 shows the map of the considered WCE lab with an 40-element SGA array (M = 8, L = 5) and vertical elements spacing BL = 3λ. We perform our numerical analysis for the point at coordinates (2,0.4,0)m on the corner of the patient walking area.
A. CASE OF ONE TRANSMITTER (K = 1) FIGURE 16 depicts the performance of GDPD, CDPD, and DOA-based localization methods. It can be easily seen that GDPD outperforms CDPD and DOA-based method by about 12 dB and 15 dB, respectively. It is obvious that for SNR = 10dB as a practical condition, the normalized location RMSE is about −13 dB (1.05 cm) for CDPD, and about −25 dB (0.07 mm) for GDPD, which confirms the drastically higher performance brought by GDPD. It can be seen that when K = 1, there is no need to use DSAR since the results are enough close to their LBs.

B. CASE OF TWO TRANSMITTERS (K = 2)
In this case, we suppose two signal transmitters, namely a WCE inside the patient's intestinal path, and a beacon  mounted on his shoulder, transmitting uncorrelated signals simultaneously. We assume about 30 cm ( 1.4λ) of vertical distance between these two transmitters, which is a harsh condition for localization. The simulation results are provided in FIGURE 17 for a WCE located at (2,0.4,0)m and a beacon located at (2,0.4,0.3)m.
It can be seen that DOA-based localization method fails while facing with more than one signal emitter. For CDPD, there is no significant changes in the results comparing to the case of having only one signal emitter. Regarding GDPD, in contrast to the case of K = 1, we observe that without applying DSAR, the GDPD is trapped to a fixed normalized RMSE of about -23 dB. This phenomenon is related to spurious peaks generated in SPF which prevent the algorithm to find the true peak on the target position. Applying DSAR technique combats the undesirable spurious peaks around the target, and brings the RMSE to its predictable values.

C. ANALYZING THE FADING EFFECT
Here, we consider that the transmitter is located at point (2,0.4,0)m and our results are averaged over 100 realizations where σ f = 0.1. FIGURE 18 shows the performance of our proposed averaging techniques to decrease the fading effect for DPD algorithm. It can be seen that VOLUME 9, 2021  the proposed averaging techniques reduce considerably the RMSE. Although the proposed methods perform almost similarly, it is worth mentioning that averaging over covariance matrix R has the lowest computational burden. Tables 1  and 2 show the complexity of operations for fading removal methods.

D. ANALYZING THE REFRACTION EFFECT
In this part, we numerically investigate the effect of refraction phenomenon on the performance of CDPD for three cases: without cancellation, pre-cancellation, and post-cancellation of refraction. FIGURE 19 depicts the value of localization RMSE normalized to wavelength (λ 21 cm), for the simulation of mentioned scenario. It can be seen that the refraction pre-cancellation method proposed in VII-B, performs such as the CDPD when there is no refraction effect. Furthermore, the refraction post-cancellation (low computational load) method proposed in VII-D, underperforms the previous technique for about 3 dB for SNR < 15 dB. For SNR ≥ 15 dB, this method cannot provide more accuracy because of the approximation error of bias cancellation. Of course, the value of −14 dB for normalized RMSE (about 1cm) is an acceptable accuracy to localize WCE for many cases. Finally, the localization method without refraction canceling is presented, whose performance is lower than two other methods as expected, where the normalized RMSE reaches −10 dB (about 2.1 cm) for SNR > 10 dB.

IX. CONCLUSION
The problem of WCE localization in WMTS labs using the combination of CDPD as a coarse (medium resolution) and GDPD as a fine (high resolution) position estimator was addressed. We showed that GDPD outperforms classical WCE localization methods in terms of accuracy (about 15dB) and capability to resolve more than one signal sources at the same time. We showed that our previously introduced technique DSAR [22] can significantly increase the robustness of GPDP while localizing several signal emitters (theoretically up to one less than the number of sensor elements in a group [13]).
We derived the CRLB formulation for SGA, and also introduced the MGLB for MGA structures, and showed that GDPD outperforms the CDPD for about 10-15 dB, depending on the location of WCE through the WMTS ward. In addition, we introduced two visual tools based on the CRLB value and also the mutual error distributions to achieve a proper geometry and elements arrangement for any arbitrary array.
Moreover, we provided a comprehensive analysis on wave refraction phenomenon in a WMTS lab due to tissue of human body, and pointed out the superiority of DPD over DOA-based localization methods in this sense by using the refraction pre-cancellation technique. We also proposed the refraction post-cancellation as a low computational load method, which of course has a lower accuracy than the former technique.
Finally, we proposed three techniques for different stages of DPD algorithm in order to combat the effect of channel fading. We particularly showed that averaging over the covariance matrix of array observations is the fastest and most effective way to reduce the destructive effects of the channel fading.

APPENDIX A THEORETICAL RMSE LOWER BOUNDS CORRESPONDED TO SGA AND MGA STRUCTURES
We derive the CRLB for DPD problem considering an arbitrary point p in a 3D monitoring area, knowing the geometry of the receiving array and the location and signal power of some determined transmitters. To derive the CRLB expression, we consider the NB matrix form of SGA structure (GDPD) in (8) which leads to the covariance matrix R in (17). Then to obtain the MGLB for MGA structure (CDPD), we replace the R with its sparse version in (19) throughout the following relations. We assume a zero mean complex Gaussian distribution for each signal snapshot, conditioned to a known parameter vector ψ. The elements of the Fisher information matrix (FIM) are given by [13] where N is the number of observations of array signals, and ψ i and ψ j are the i-th and j-th elements of parameter vector ψ, respectively. To comply with 3D applications like WMTS and by assuming that fading coefficients are constant during one observation interval, we consider the unknown parameter ψ = p = [x, y, z] T . Assuming that noises at receiver elements are uncorrelated and have the same power σ 2 , we can rewrite (17) as from which the elements of F can be determined. To calculate R −1 , according to the matrix inversion lemma we have For a single transmitter (K = 1) with power P, from (50) we get On the other hand, derivation of (49) with respect to ψ ∈ {x, y, z} leads to To compute the CRLB for localization of the k-th transmitter among K ones, i.e., p = p k , we notice that all transmitters except the one at point p k are considered to have stationary locations and so, we can ignore them during the derivation procedure. Then we have ∂ã k ∂ψ , . . . , 0 . (53) For notation simplicity, we remove the subscript k to havẽ a =ã k . Assuming that the emitted signal from p is uncorrelated to other emitter signals, we have where P is the k-th diagonal element of R S representing the transmitted power from point p that we assume it is normalized to one.
Considering (3), for the u-th element ofã we havẽ where d u is the distance between p and the u-th receiver element with where φ ψu is the angle between ψ-axis and the line pointing p from the u-th element. We have where ζ = [ζ 1 , . . . , ζ U ] T , ζ u = −(1/d u + j2π/λ), φ ψ = [φ ψ1 . . . , φ ψU ] T . Using (54) and by analogy, we get where C ψ diag(ζ • cos φ ψ ) 1 U + 1 U diag(ζ * • cos φ ψ ). (60) For a CDPD problem with M number of L-element independent arrays, we have where W = I M ⊗ 1 L . In fact, W masks any statistical dependency between observations at different arrays, which leads to a tighter LB due to the lack of mutual information between different arrays. Finally, we can compute the elements of F and consequently derive the CRLB or MGLB as tr(F −1 ). In this regard, we notify some important points, as follows.
• For an unbiased position estimatep, we have where rms is the RMS of localization error.
• Due to the fading coefficients included inã, which are considered constant during one DPD estimation interval, we can average the LBs resulted from several random sets of fading coefficients to achieve a more reliable LB [29]. VOLUME 9, 2021 • The total SNR in (49) is calculated as Moreover, the specific SNR related to the k-th transmitter, i.e., SNR k , by using the properties of tr() operator, as is given where P k is the power of the k-th signal transmitter.

APPENDIX B LOCALIZATION BIAS DUE TO REFRACTION
Considering the scenario depicted in FIGURE 11, we investigate the virtual displacement of the arbitrary point p = p u inside a circle with radius R due to the refraction phenomenon. The wave radiated from p hits the perimeter of the circle at point r with angle θ and leaves r with refracted angle θ according to the Snell's law, with refraction indices are n > 1 inside, and n = 1 outside the circle, respectively. Therefore the point p with distance l from the refraction point r will be seen at the point p with distance l from r. Now we find the displacement value ρ as a function of ϕ which is the central angle of refraction point r. We have Using some geometric rules, it can be easily shown that Considering (65) and specially the definitions in (66), and using some trigonometrical manipulations, (67) can be written as FIGURE 20 depicts ρ(ϕ) for three different values of d, which is the distance between p and the center of circle. In the first case (rightmost) where d = 0, we get a circle that is concentric with body circle, which makes no localization bias. In the second case, where d = 0.3R, we get a convex connected shape with a little bias or distance between real location and estimated one as shape's centroid. The third case, where d = 0.5R, generates a concave closed shape with larger bias compared to the first one. Finally, the last case with d = 0.8R produces a concave open curve with a very larger bias than the two former cases. The open parts of curve are related to the situations where the emitter is located at hidden areas as described in Section VII-C and depicted in FIGURE 13.
Finally, assuming the array symmetry such as in FIG-URE 11, the refraction intrinsic bias in (39) equalsd − d, whered is the distance of the centroid of p locus, from the cylinder axis. Unfortunately the integral of the function in (68) cannot be derived analytically with respect to ϕ, so we should use the numerical averaging to calculate the estimation bias.