**2. Overview of Approach**

Techniques based on cross-correlation and generalized correlation (GCC) [8] have been employed to determine the time difference of arrival of the signals (TDOA) given its computational cost and accuracy of the results. To obtain a better estimate of the TDOA, it is convenient to filter the signal before its integration, as shown in Figure 1 [6].

**Figure 1.** Scheme for obtaining the time of arrival (TOA).

The cross-correlation *Rxixj* between the signal *xi* and *xj* filtered by the filters *Hi* and *Hj*, is expressed as a function of the power spectral density *Gxixj* , as shown below.

$$R\_{\mathbf{x}\_i^x \mathbf{r}\_j^\circ}^{\rm GCC}(t') = \int\_{-\infty}^{+\infty} H\_i(f) H\_j^\circ(f) \mathbf{G}\_{\mathbf{x}\_i \mathbf{r}\_j^\circ}(f) e^{j2\pi ft'} df = \int\_{-\infty}^{+\infty} \varphi^{\rm GCC}(f) \mathbf{G}\_{\mathbf{x}\_i \mathbf{r}\_j^\circ}(f) e^{j2\pi ft'} df \tag{1}$$

where {\*} indicates a conjugated complex and ϕ*GCC*(*f*) is a frequency-dependent weight function. Due to finite observations, we can only obtain an estimation of *Gxixj* (*f*) [9]. Therefore, to obtain the TDOA, the following expression will be used [6].

$$\mathcal{R}\_{\mathbf{x};\mathbf{x}\_{j}}^{\rm GCC1}(t') = \int\_{-\infty}^{+\infty} q^{\rm GCC}(f) \mathcal{G}\_{\mathbf{x};\mathbf{x}\_{j}}(f) e^{j2\pi ft'} df \tag{2}$$

where *<sup>G</sup>*<sup>ˆ</sup> *xixj* (*f*) is the obtained estimation of *Gxixj* (*f*). For each pair of sensors, the TDOA is taken as the time delay that maximizes the cross-correlation between the filtered signals of both sensors, that is: τˆ*GCC ij* <sup>=</sup> arg max*t*- *<sup>R</sup>*<sup>ˆ</sup> *GCC*<sup>1</sup> *xixj* (*<sup>t</sup>* - ) .

A general model for three-dimensional (3-D) estimation of a source using *M* receivers is developed. To obtain the location of the source, we start by knowing the spatial position (*xi*, *yi*, *zi*) of a certain number of sensors. Let (*xs*, *ys*, *zs*,), the position of the source to be located, the distance between the source and the *i*-th sensor will be:

$$d\_i = \sqrt{(\mathbf{x}\_i - \mathbf{x}\_s)^2 + (y\_i - \mathbf{x}\_s)^2 + (z\_i - \mathbf{x}\_s)^2} \tag{3}$$

The range difference in distance between the *i*-th receiver and the first receiver, *di*1, is given by:

$$d\_{l1} = c \cdot \tau\_{l1} = d\_l - d\_1 = \sqrt{(\mathbf{x}\_i - \mathbf{x}\_s)^2 + (y\_i - y\_s)^2 + (z\_i - z\_s)^2} - \sqrt{(\mathbf{x}\_1 - \mathbf{x}\_s)^2 + (y\_1 - y\_s)^2 + (z\_1 - z\_s)^2} \tag{4}$$

where *c* is the sound velocity in the medium, *d*<sup>1</sup> is the distance between the first receiver and the source, and τ*i*<sup>1</sup> is the estimated TDOA between the first receiver and the *i*-th receiver [10].

Equation (4) considered for all the sensors form a nonlinear equation system whose solution can be found by several ways. After studying different resolution methods, it was decided to use the Newton-Raphson method since it offers very good results and computation time.

### *Newton-Raphson Method*

To get the position of the thermoacoustic source inside the medium, we have solved the nonlinear system using the Newton-Raphson method [11] by means of partial derivatives. Consider a system of *m* equations and *n* unknowns.

$$f\_{\mathfrak{m}}(\mathbf{x}\_1, \mathbf{x}\_2, \mathbf{x}\_3, \dots, \mathbf{x}\_{\mathfrak{n}}) = 0 \tag{5}$$

This system can be written in vector form as *f*(*x*) = 0, where *f* is a vector of *m* dimensions and *x* is a vector of *n* dimensions. To solve this system of equations, we have to find a vector *x* such that the function *f*(*x*) equals the null vector. If we call η to the solution of the system and *xr* to an approximation of it, we can develop *f* in Taylor series around this approximation as:

$$f(\mathbf{x}) = f(\mathbf{x}\_r) + \nabla f(\mathbf{x}\_r) \cdot (\mathbf{x} - \mathbf{x}\_r) + \dotsb \tag{6}$$

Since *f*(η) = 0 then we get, as an approximation, the following.

$$0 \approx f(\mathbf{x}\_r) + \nabla f(\mathbf{x}\_r) \cdot (\eta - \mathbf{x}\_r) \tag{7}$$

Now, we can define the vector *xr*+<sup>1</sup> as this approximation, which is closer to the root than *xr*. We can continue with the iterative method to obtain approximations closer and closer to the solution. To write the iterative method, the term ∇*f*(*xr*) is replaced by the Jacobian of the function *f*, that is:

$$f(\mathbf{x}\_{\tau}) = \begin{bmatrix} \frac{\partial f\_1}{\partial x\_1} & \frac{\partial f\_1}{\partial x\_2} & \dots & \frac{\partial f\_1}{\partial x\_n} \\ \frac{\partial f\_2}{\partial x\_1} & \frac{\partial f\_2}{\partial x\_2} & \dots & \frac{\partial f\_2}{\partial x\_n} \\ \vdots & \vdots & \vdots & \vdots \\ \frac{\partial f\_m}{\partial x\_1} & \frac{\partial f\_m}{\partial x\_2} & \dots & \frac{\partial f\_m}{\partial x\_n} \end{bmatrix} \tag{8}$$

which is a *n* × *m* matrix. Then, it is possible to obtain a new value of *xr*+<sup>1</sup> by solving the following relationship.

$$f(\mathbf{x}\_r)(\mathbf{x}\_{r+1} - \mathbf{x}\_r) = f(\mathbf{x}\_r) \tag{9}$$

Then, iteratively, we can approximate more and more the *xr*+<sup>1</sup> to η until a solution error *xr* <sup>−</sup> *xr*+<sup>1</sup> previously fixed is reached.

This method has been compared with other algorithms for solving systems of nonlinear equations and has shown some advantages over the rest like good accuracy and low computing cost. However, if the initial solution value of the system differs greatly from the real solution, then the method does not converge conveniently. Figure 2 shows the convergence of the location algorithm as a function of the distances between the initial point of the method, the result of the reconstruction of the position and the real position for a total of 10,000 simulations.

*Sensors* **2019**, *19*, 1971

**Figure 2.** (**a**) The distance between the initial position for the algorithm and the real position of the source is shown on the abscissa axis, while the axis of the ordinates shows the distance between the signal reconstructed by the algorithm and the real position of the source; (**b**) The reconstructed positions for each of the 10,000 simulations are shown together with the real position (black).
