Distance Estimation Algorithm for Wireless Localization Systems Based on Lyapunov Sensitivity Theory

The paper describes a novel distance evaluation algorithm based on the time-difference of arrival (TDOA) principle. The proposed method solves the distance estimation problem applying the Lyapunov theory. To perform this task, the distance evaluation problem is converted to a parameters identification process exploiting the sensing signal peculiarities. This latter combines the properties of the Frank-Zadoff-Chu (FZC) sequences with the Orthogonal Frequency-Division Multiplexing (OFDM) modulation scheme. The resulting signal can be used to improve accuracy and precision in distance estimation using a reduced signal bandwidth. The system was modeled considering either an additive white noise Gaussian channel and a moderate multipath model. The simulation results demonstrate that the system achieves an accuracy and a precision better than 10 cm considering SNR equal to 0 dB and signal bandwidths of only 125 MHz.


I. INTRODUCTION
Wireless localization systems are extensively used in modern communication applications to furnish different services or to execute several tasks. The main issue for these systems is to achieve the precision and the accuracy required by the target application. For this reason, it is very hard to implement a general solution for such a variety of applications [1], [2].
An important classification for wireless positioning systems can be made considering the measurement principle employed to evaluate the target distance and, thus, the position. Three different measurement methods are used: angle-of-arrival (AOA) or phase-of-arrival (POA) [3], [4], received signal strength (RSS) [5], [6], and propagation-time based system that can further be divided into three subsets: round trip-time-of-flight (RTOF), time-of-arrival (TOA) and time-difference of arrival (TDOA) [7].
Focusing the attention on the last category, the TDOA-based systems are able to evaluate the target position by computing the time-difference of arrival of at least two received signals. To this purpose, the object under control is assumed to be a transmitting unit and at least two receiving The associate editor coordinating the review of this manuscript and approving it for publication was Fang Yang . units are needed. In these assumptions, the TDOA system requires only a synchronization between the measuring units (receivers) but not with the target unit (transmitter). The physical layout calls for a backbone network or a reference transponder in a known position, in order to correctly calibrate the system. Basically, the measuring units, in the propagation-time based systems, are able to evaluate only the distances from the target unit. Based on several of such measurements, a 2-D or 3-D position can be derived using lateralization or triangulation methods [8], [9].
The literature proposes several approaches to evaluate the TOA/TDOA. The most accurate methods use Ultra-Wide-Band signals (UWB). For example, Wang et al. [10] propose a radar sensor that evaluates the TOA of a static or slow-moving target combining the advantages of the interferometry technology with the properties of the Frequency-Modulated Continuous Wave (FMCW) signals. The radar system operates at 5.8 GHz ISM band and exploits a signal of 160 MHz, achieving an average accuracy of 1.65 cm. Another example of FMCW system is presented by Waldmann et al. [11]. In their work the FMCW signal is combined with short pulses with a bandwidth of 1 GHz centered to 7.5 GHz to achieve an error of 1.7 cm in absence of multipath distortions. Mahfouz et al. [14] proposes a high-accuracy indoor 3D positioning system that use an 8 GHz carrier signal modulated by a 300-ps Gaussian pulse. The Mahfouz system exploits a GPS-like scheme with four receivers that localize the position of a mobile transmitter. The estimation process involves the detection of the first arrival path of the pulse and it allows to distinguish the paths related to the multipath delays. The impulse-radio is generated in time domain and is shaped and filtered in order to obtain the wanted frequency spectrum. The experimental results demonstrated that the system achieves a precision of 1.53 mm in 3D case, with an accuracy of only 3.62 mm and a worst case of 5.6 mm. Lipka et al. [17] proposes a novel localization concept whereby the phase of a 24GHz-FMCW signal with a bandwidth of 250 MHz is fed into an Extended Kalman Filter without any preprocessing. Their experiments in severe multipath conditions show that the system is able to locate the target with a 3D accuracy of 1.7 cm.
Another technique to measure the distance involves the transmission of a Multi-Carrier UWB signal (MC-UWB) [18]. In these systems, the transmitted signal comprises many subcarriers, orthogonally spaced. They don't exhibit the time-domain resolution offered by the Impulse-Radio UWB systems (IR-UWB), but allow to evaluate the phase of arrival to enhance the system precision and accuracy. The MC-UWB systems occupy a narrower bandwidth in comparison to the IR-UWB systems and, most of all, they offer a high spectral efficiency and flexibility. In fact, the subcarriers can be nullified or placed at a frequency that allows to co-exist with other systems occupying the same frequency-band. Ayhan et al. [15] and Redant and Dehaene [16] proposes a sub-mm ranging system, which estimates the TOA of a RF signal between two nodes with a frequency domain estimator. In their work, a simulation precision better than 1 mm is obtained for SNRs below 0 dB, by transmitting an Orthogonal Frequency Division Multiplexing (OFDM)-like signal whose duration is 9 µs, with a 6 GHz bandwidth on a 60 GHz carrier. The ranging distance is divided in three main steps: in the first one a coarse distance value is obtained evaluating the correct starting point of the receive signal. In the second one the previous value is fine adjusted through a frequency domain estimator that evaluate the phase of the received signal. In the last one a super fine range estimation is performed calculating the carrier phase offset introduced by the TOA of the system.
Another important advantage of the MC-UWB systems is that the baseband signal in the transmitter path can be generated in digital domain and then with a Digital to Analog Converter (DAC), converted in analog domain. In this manner, the transmitter architecture is very simple, as the orthogonality between the subcarriers can be obtained applying the Inverse Discrete Fourier Transform operation in digital domain. Moreover, in this domain any offset in the signal amplitude, phase or time-delay can be easily avoided or recovered.
As mentioned, the use of an UWB signal allows to extract the target distance with a very high accuracy and precision. However, the system design is more complex because the performance of the analog blocks must be guaranteed for huge signal bandwidths. Moreover, the digital part should process a large amount of data and the overall cost increases accordingly.
The aim of this paper is to introduces a novel TDOA algorithm that is able to estimate the distance with a very high accuracy and precision, using small signal bandwidths. Compared with the most accurate methods in literature, the proposed approach allows achieving the same distance error but with a narrower bandwidth. The resulting system, thus, is less expensive than the others approaches that exploit Ultra-Wide signal bandwidth.
The method is based on the Lyapunov theory and transforms the distance evaluation problem in a parameter identification process. The TDOA is derived in frequency domain from the channel transfer function, H (f ), related to the difference between the received signals. This latter is obtained exploiting the properties of the sensing signal that combines the Frank-Zadoff-Chu (FZC) mathematical sequences with the OFDM modulation scheme [19]- [24].
The parameter identification process is brought out by evaluating the square of the prediction error. The idea is to design a stable artificial dynamic model characterized by an equilibrium point coinciding with the minimum of the error prediction [22].
The paper is organized as follows: section II introduces the channel model used to derive the programming problem and describes, in details, the mathematical formulation of the distance evaluation algorithm. Section III presents the simulation results derived considering in the first set of tests a ideal channel composed by a single ray. Then a multipath channel is considered to evaluate the algorithm performance in terms of precision and accuracy even for moderate harsh environments. The conclusions will close this work.

II. ALGORITHM PRINCIPLES A. SYSTEM MODEL
Let's define the transmitted signal, x(t). In our system x(t) represents an OFDM symbol composed only by pilot subcarriers. Each of them representing a coefficient of a Frank-Zadoff-Chu (FZC) mathematical sequence [19], [20], [24]. The ZC sequences are part of the Constant Amplitude Zero-Auto-correlation family and, due to their properties, they are massively used in telecommunications systems to resolve synchronization problems [25], [26].
Assuming h(t) as the channel impulse response, the baseband received signal is obtained from the convolution between x(t) and h(t): where n(t) is zero-mean additive white Gaussian noise introduced by the channel. Due to the properties of the Fourier transform, the spectrum Y (f ) can be obtained by the multiplying the transmitted signal spectrum, X (f ), times the channel VOLUME 7, 2019 transfer function, H (f ): where N (f ) = F {n(t)}. If we multiply the eq. 2 by the complex conjugate of X (f ), we obtain the expression below: The transmitted signal and the noise process are statistically independent. The second term, in the expression above, can be neglected. On the contrary, in the first term the square modulus of X (f ) is equal to 1 as the signal has a constant amplitude. It's clear that the result of the multiplication is equal to the channel transfer function, H (f ). No assumptions about the signal origin was made and, hence, the received signal could either be directly emitted by a wireless transmitter (Line-Of-Sight, LOS) or be a component reflected by any object in the operating scenario (Non-Line Of Sight, NLOS). Generally, the channel impulse response can be modeled as [27], [28]: where χ(t) represents the distortion of the i-th arriving component due to the frequency selectivity of the interactions with the environment, α i is the amplitude related to the i th signal component and τ i the time delay associated with the i th reflection and δ is the Kronecker function. N is the number of scattered components, and * represents the mathematical operation of convolution. The distortion functions are in general unknown and, thus, a simplified channel model can be used: with L > N . The time delay τ i depends by the distance d i between the transmitter and the receiver (or between the receiver and the specific obstacle) and the signal propagation velocity, v, (in air is around 3 × 10 8 m/s): A distorted pulse looks like a sequence of closely-spaced pulses with amplitude determined by the sum of the power carried by the multipath components as well as the pulse distortion. Thus, the simplified model might identify a number of scatters in excess to those physically existing, thus generating the so called ghost components but their locations will be closely spaced around the locations of the true scatters.

B. MATHEMATICAL FORMULATION
The eq. 5 can be obtained by the inverse Fourier transform of eq. (3). Under this statement, the implemented signal has a pulse correlation function. This result can be exploited to evaluate the value of the time delays τ i , that represents the system TDOA. It's important to note that in the model expression the Carrier Frequency Offset (CFO) was neglected and, thus, we will consider the baseband signal where the phase and the frequency of the local oscillation are perfectly recovered.
Let's consider the channel model in (5), if we apply the Fourier transform, the resulting transfer function is the sum of L complex sinusoidal waves, with a proper value of amplitude α i and time delay τ i : The transfer function H (f ) can be adopted as estimate model to evaluate the parameters α i and τ i and, from these, the distance between the transmitter and the receivers. To obtain this objective a parameter identification process can be implemented. The parameters are extracted from the available data concerning a given system using the prediction error concept. The prediction error is the difference between a predicted output signal and the measured input signal. The former includes an estimation model, that measures the reference model and relates the measured output to its unknown parameters. In our case, the complex nature of the transfer function expressed in (7) implies a prediction error expressed by two components, real and imaginary as follows: where Y r and (Y i ) are the reference models coincident to the model dataset. The latter is obtained from the received signal from which we aim to extract the parameter information associated with the estimated model. The functions y r (f , α i , τ i ) and y i (f , α i , τ i ) are the predicted signal including the unknown parameters, expressed as follows: The aim of the parameter identification process is to evaluate the parameters included in estimated model which nullify the prediction error vector: The proposed procedure aims to solve the problem in (11), minimizing the following scalar positive semi-definite function: The traditional approaches try to solve the first-order derivative expression: However, the problem can be not easily solved due to the non-linear nature of the set of equations. In fact the traditional numeric methods, for example the Newton-Raphson, require matrix inversions and/or factorization.
The suggested procedure overcome this problem. The method is based on the idea presented in [21] and [22] and it comprises the following three principles: • conversion of the generic multi-objective mathematical programming problem in a set of differential equations; • scalarization of the multi-objective optimization problem by a definite positive quadratic form; • generation of a computational paradigm based on Lyapunov theory [29]; The set of differential equation must be simultaneously satisfied and the last point brings the design of a stable artificial dynamic model characterized by an equilibrium point equivalent to the minimum of W (f , z).
The components of the vector z are assumed to evolve according to a proper time dynamic, z i (t * ). This assumption leads to consider the eq. (12) as a Lyapunov function. Consequently, if its time derivative is forced to be negative-definite or negative semi-definite negative along the trajectory z i (t * ), then the Lyapunov theorem dictates the existence of an asymptotic stable equilibrium point minimizing W (f , z) and, consequently, a solution for the parameter identification process [29]. In order to prove this statement let's compute the time derivative of W (f , z) applying the chain rule to the eq. (12): The time derivative of the composite error function can written asĖ whereż is the time-derivative of the parameters. Then, by substituting the equation above in (14), W (f , z) can be expressed as follow:Ẇ Therefore, if z(t * ) evolves according to the gradient of W (f , z): where k is a positive constant that may assume a different value for the two state variables (in this case k = |k α , k τ | T ). Then by combining the eq. (17) with (15) and (16), we obtain: (19) Note that the latter equation has a quadratic form and it is certainly negative-semidefinite. Due to this important result, it is possible to postulate that if z(t * ) evolves according to the eq. (17), then the Lyapunov conditions are satisfied. In facts, W (f , z) is positive-semidefinite and the time derivative is negative-semidefinite. Under these conditions, the asymptotic stability of the equilibrium point is guaranteed. Moreover, by observing that the matrix (dW (f , z)/dz) (dW (f , z)/dz) T is symmetric and positive-definite and its eigenvalues are real and positive, it follows that E(f , z) converges to the equilibrium point E(f , z * ) = 0 exponentially.
The proposed scalarization of the problem guarantees that the solution obtained is a feasible Pareto solution [23].
In this dynamic system, the variable t * is an artificial parameter that allows to observe the transient nature of the trajectory z. Hence, there is no relationship between t * and the real time instant, t, of the received/transmitted signal. More in details, the eq. (19) can be re-written as follows: The solution of the artificial system in (17) is equivalent to the solution of the problem in (13). It is important to highlight that such solution doesn't require any matrix inversion and/or factorization, overcoming the main difficulties arising by the ill-conditioning of the Jacobian matrix. This assumption represents an important advantage of the proposed method because it ensures the presence of a problem solution.
C. ALGORITHM DESCRIPTION Figure 1 sketches the block diagram of the proposed distance estimation algorithm, divided in two main steps. In the first one the algorithm presented in [30] is used to calculate the initial state vector z 0 . The latter affects the transient dynamics of the artificial system in (20). Obviously, an initial state close to the equilibrium point increases the speed of the convergence. Moreover, due to the nature of the distance evaluation problem, the components number, L, in the channel model is not a priori-known because it depends by the operating environment. The vector z 0 can be evaluated from the cross-correlation function between the received signal, y(t), and its transmitted local copy, x(t). To perform this task, the Fourier Transform is applied on y(t) and x(t) to obtain the related signal spectrum, Y (f ) and X (f ) respectively. Then, the conjugate copy of the latter, X * (f ) is multiplied by Y (f ) [30]. This term represent the dataset required by the parameter identification process. Applying the Inverse Fourier Transform on Y (f )X * (f ), the cross-correlation function can be obtained. As mentioned, it is a pulse function in which the pulses, of amplitude α i , are located at the time τ i .

FIGURE 1. Distance evaluation algorithm -block diagram.
A search algorithm is implemented to evaluate these peaks and the result is the state vector z 0 . It's important to note that the correlation-based algorithm is able to evaluate only the integer part of the i-th time delay τ i that will be a multiple of the sampling time. Then, from the simulation of the artificial dynamic system it is possible to adjust the fractional part of the delay, τ i . The algorithm 1 resumes the procedure described above.

Algorithm 1 Estimation Algorithm for Finding the Start
Conditions of the Artificial Dynamic System 1: function StartConditionsestimation (y(t), X (f )) 2: Input: y(t) → received signal X (f ) → transmitted signal spectrum 3: Output: z 0 → path gain and time-delay vector 4: In the second step, the initial state vector z 0 , and the reference model, H (f , α, τ ), are passed to the algorithm to compute the prediction error and simulate the artificial dynamic system in (20). Once the artificial dynamic system reaches the equilibrium point, the algorithm returns the state vector z containing the estimated values α and τ . The algorithm in 2 reports the pseudo code of the procedure, whereas in fig. 2 is reported the block diagram of the artificial system for the state variable z. An important feature characterizing the proposed paradigm is the intrinsic data filtering capability derived by the integral action of the computational process and by the use of the cross-correlation function used to estimate the initial state vector. This makes the proposed approach particularly powerful in solving the optimization problems in presence of uncertainty or noisy data.
The gain factor, k, plays an important role in the convergence of the method. Small values of k usually lead to an inaccurate solution. Generally speaking, the accuracy of the solution will be improved by increasing k within a small range. But beyond some point, any further increase of the gain factor may result in slower convergence, and possibly, oscillations. This phenomenon is similarly to what happens in the gradient search method in optimization. The trend of the solution as a function of the gain, k, should be used during the simulation process in order to empirically choose their optimal values. An important advantage of the proposed method is the possibility to include a saturation block within the artificial dynamic system (fig. 2). This feature allows to include constraint relationships related to any physic problem, within the simulation of the artificial dynamic system. Let's consider the problem of remote-distance estimation in outdoor environment with the presence of two active devices: a transmitter and a receiver pair. In this scenario, a LOS path exists between the active devices and it is possible to assume that the values of τ i decrease monotonically because the delays grow with the distance. A similar assumption can be made for the amplitude α i when a LOS system is considered. As a matter of facts, the presence of an active transmitter in the system ensures an inverse quadratic distance dependence for the echo amplitude [31]. Starting from these assumptions we can re-formalize the problem (11) to include the constrains equations:

Algorithm 2 Parameters Identification Algorithm
The proposed procedure can be applied considering this new problem formulation.

A. IDEAL CHANNEL
To validate the proposed distance evaluation algorithm, several simulations were performed considering different values for the bandwidth and the sampling time.
In the first simulation set, an ideal channel composed by a single ray was considered (L = 1). The echo represents the direct path between the system nodes. The channel is assumed to superpose an additive white Gaussian noise to the signal and the performance were extracted as a function of the SNR at the receiver side. The signal comprises 128 sub-carriers allocated in a bandwidth of 125 MHz and sampled with a 250 MHz sampling frequency to obtain a total of 256 representation points. It was modeled a transmitter movement along the line connecting the first transmitter position with the receiver, considering a distance ranging from 1m up to 18m. In this test-bench we focus our interest only on the channel effects so we can assume, for sake of simplicity, that the receiver and the transmitter antennas have the same directivity patterns and the phase center is independent within the used bandwidth, and that other impairment sources related to the circuit hardware can be neglected. Due to the different values for the amplitude vector α and the time-delay vector τ two different gains k α and k τ were considered (k α = 0.001 and k τ = 100). The differential equations describing the artificial dynamic system was solved using the ode15s solver implemented in MATLAB.
In fig. 3 is represented an example of the trajectories of α(t * ) (normalized value) and τ (t * ) (measured in samples) as function of time t * for both the error functions E 1 and E 2 . Note that after a transient time of 18 ns the artificial dynamic system associated to τ (t) reaches an equilibrium point whereas for α(t * ), the system equilibrium is achieved after 6 ns.
As previously mentioned, the transient time depends by the k τ gain value. In fact, in fig. 4 it is reported an example of the    E (f , z).
influence of k τ on the equilibrium state for the considered dynamic system. For a value of 50, the system reaches its equilibrium after 15 ns whereas with a value of 100 the transient ends in 6 ns. Note from figure 4 that for a k τ = 200 the system achieves its equilibrium state within a time greater than that in the k τ = 100 case. This fact is due to the larger oscillations introduced by a greater value of k τ .
The fig. 5 sketches the trends of the square modulus of the error function as a function of the artificial parameter t * . Once the equilibrium was reached, the problem in (21) was solved achieving a final value of 8.60e-8 for the function E 1 and 5.7e-5 for E 2 . This discrepancy in error values is due to the different gains and simulation times considered for the functions.
In fig. 6 are reported a set of estimated distances for a SNR of 0 dB and SNR of 15 dB compared with the values extracted from the cross-correlation function R xx (t). As previously mentioned, the accuracy and the precision of distances extracted from R xx (t) are affected by the sampling time. To enhance the precision and the accuracy in correlation-based methods, the sampling time must be accordingly increased. This assumption represents a major drawback of the correlation-based methods. In fact, increasing the sampling time, the system complexity grows as well. Moreover, the hardware must process a much greater amount of data exhibiting higher costs when compared to the circuits needed for other methods that are not based on the correlation function.
The proposed solution utilizes the cross-correlation function to compute only the state vectors for the artificial dynamic system. If the Nyquist theorem is satisfied, the artificial dynamic system is able to evaluate the complete value of the time-delay correctly, as shown by the linear curve in fig. 6. Hence, the proposed solution is less sensitive to the sampling time reducing the costs and the performance required by the hardware.
In fig. 7 are reported the performance of the proposed distance evaluation algorithm in terms of accuracy and precision and as a function of the SNR. The performance were evaluated considering a signal composed by 128 sub-carriers allocated in 125 MHz of signal bandwidth. The signal is sampled with a sampling time of 250 MHz. The system accuracy was compared with the Cramer Rao Lower Bound (CRLB). For a single-path AWGN channel, the CRLB is given by [32]: where v p is the propagation velocity (3e8 m/s in air), and β is the effective signal bandwidth. The CRLB for time-based ranging decreases with the square-root of the SNR and effective signal bandwidth. The graphs evidence that the proposed solution is able to extract the distance with a mm precision and accuracy for a SNR greater than 15dB. For a SNR lower than 15 dB, obviously, the performance are reduced but the error is lower  than 10 cm even if the received signal has a 0dB SNR. This is a rather extreme condition in a common scenario. In fact, the proposed system is based on an active transmitting target and even for a low radiated power of the order of −10dBm @ 5GHz and a distance of a hundred of meters ensures a LOS received signal level in excess to −60dBm.
In fig. 8 are reported the system performance as a function of the signal bandwidth and of the number of sub-carriers. From the figure, it may be noted that increasing the bandwidth or the number of sub-carriers, the system precision and accuracy increase as well. Doubling the bandwidth is equivalent to double the number of sub-carriers in terms of system performance. In facts, if the system transmits a signal with 256 sub-carriers allocated in 125MHz of bandwidth, it reaches the same precision and accuracy of a system that uses a signal composed by 128 sub-carriers allocated in a 250 MHz bandwidth. This linear increment is another advantage of the proposed method, as the values of these two parameters can be selected to achieve the best trade-off between cost, hardware complexity, convergence speed of the algorithm and so on.

B. MULTIPATH CHANNEL
In channels affected by multipath, the proposed algorithm must evaluate the amplitude and the time delay associated to the single sinusoidal wave representing the specific signal echo. In the first test-bench a channel composed by two echoes was considered. The first one represents the direct path between the system nodes and the second one is a multipath component. For the direct echo a fixed distance d 0 of 5 meters was imposed, while, a variable distance in the d 0 + 0.9m and d 0 + 5mm range was considered. Due to signal bandwidth and the system sampling time, echoes placed beneath to the breakdown distance of d 0 + 0.9 are unresolved by the system because the cross-correlation function will exhibit a single pulse that contains both echoes [35]. Further, a LOS scenario was simulated, so the direct path is the first arriving echo at the receiver side, and, it exhibit the strongest amplitude (normalized to 1), whereas, a normalized amplitude of 0.8 was chosen for the multipath echo to simulate the arrival of a very strong multipath interference.
In fig. 9(a) are represented the components of α(t * ) and τ (t * ) evaluated by the artificial dynamic system,  setting k α = 0.002 and k τ = 200. Note that for both components the system reaches the equilibrium point after 6 ns. In fig. 10 an example of the error functions for both E 1 and E 2 are represented in the case of multipath channel. Once the equilibrium was reached, the problem in (21) was solved achieving a final value of 9.7e-8 for E 1 and 7e-4 for E 2 . Fig. 11 shows the system performance in terms of precision and accuracy for both the direct signal and the multipath component. The performance are evaluated transmitting a signal composed by 128 sub-carriers arranged in a 125MHz bandwidth and sampled with 500 MHz of sampling frequency. The accuracy is compared with the CRLB that is slightly different than the single path AWGN channel case. In facts, it must be evaluated for each of the delay values τ i , using the expression below [32], [33]: where E s is the energy of the signal, x(t), given by: and E s and E s are its time derivatives. Due to the properties of the FCZ sequences, the signal energy is constant and the expression of the CRLB in (23) depends on the signal parameters and the path gains. As expected, the performance are slightly worse than in the ideal case. In fact, considering 15dB of SNR, the precision in multipath case is equal to 1.8 cm while in ideal case it was around 1 cm. Analog considerations are valid for the accuracy, in facts, ever considering a 15dB SNR, in the ideal case it was 1.1 cm, whereas, in multipath conditions the direct echo estimate is 2.3 cm. Although the performance in multipath channel evaluation are lower than in the ideal case, they represent a good level of precision and accuracy even in presence of a strong multipath interference. Hence, the algorithm is very robust against the multipath problem. As mentioned, to further enhance the system performance the bandwidth or the number of sub-carriers can be changed to achieve the target precision and accuracy as function of the multipath echoes present in the channel. Moreover, it is important to note that the direct echo was estimated with a very high level of accuracy and precision directly from the channel model in (7). Any equalizer filter or iterative methods, as for example the search and subtract approach [27], [35], was introduced in the system.
Finally, in the last set of test-benches a moderate harsh outdoor environment was simulated. The channel model comprises the LOS echo and 4 dominant multipath components and the transmitted signal comprises 128 sub-carriers allocated in a 125 MHz bandwidth and sampled with a sampling frequency of 500 MHz. To extract the components of α(t * ) and τ (t * ) a value of 0.002 for k α and 200 for k τ were imposed. In fig. 12 an example of α(t * ) and τ (t * ) components are reported. Note that the system is able to achieve a stable equilibrium point for all vectors component in 5 ns.  The trend of the error E(f , z) as a function of the artificial parameter t * is sketched in fig. 13 for both E 1 and E 2 . Once the equilibrium was reached, the problem in (21) was solved achieving a final value of 1.37e-7 for E 1 function and 1.3e-3 for E 2 .
In fig. 14 the final performance of the algorithm in presence of moderately harsh outdoor environment are reported. Increasing the multipath components' number, the accuracy and the precision decrease. Considering a SNR of 15dB the proposed algorithm is able to evaluate the distance with a precision and an accuracy better than 10 cm. However, the system achieves a good level of accuracy and precision even for lower SNRs.
In table 1 the accuracy of the proposed method is compared with the most accurate approaches to evalute the target distance present in literature. The methods reported in the table exploit Ultra-Wide Band signal (B W ≥ 500MHz) and they are able to compute the target distance in mm or sub-mm range. As reported in fig. 14 in the case of moderate harsh environments, the proposed method achieves an accuracy better than 10 cm even for a SNR of 0 dB. By comparing the proposed approach with the other methods reported, we can affirm that the system is able to achieve the same performance in distance estimation using smaller signal bandwidths (Bw = 125MHz). Moreover, the performance can be further improved implementing digital techniques to enhance the signal-to-noise ratio at the receiver side or implementing digital filtering as for example the Kalman Filter.

IV. CONCLUSION
A novel approach to evaluate the TDOA between a measuring units pair and an active target was presented. The method is based on the Lyapunov theory and allows to convert the distance estimation problem into a parameter identification process. The sensing signal is a Multicarrier-Wide-Band signal composed by the coefficients of the Frank-Zadoff-Chu mathematical sequence. Its properties are exploited to estimate the distance with a very high accuracy and precision. In facts, in moderate harsh environment the proposed procedure reports a simulated error lower than 10 cm, considering a signal composed by 128 sub-carriers allocated in a 125 MHz bandwidth. The accuracy and precision achieved are comparable with the most accurate and precise approaches in literature. However, the proposed method uses a signal with a much narrower bandwidth with respect to the other approaches that, vice versa, use UWB signals.
The system hardware is based on the Software-Defined Radio architecture and the signal parameters can be adjusted to evaluate the algorithm performance as a function of these two parameters. In a future development, the carrier frequency offset will be included in the channel model and model itself will be validated for outdoor/indoor environments.