## 3 Fast Hyperbolic Hough Transform
### Fast HHT
Thus, using (21) and (22), in the tile with center \(P\), an upper bound \(\overline{\mathcal{L}}_{n}\left(P\middle|\widetilde{d}_{i,1}\right)\) for \(\mathcal{L}_{n}\left(P^{\prime}\middle|\widetilde{d}_{i,1}\right)\) can be obtained, as follows:  
\[\overline{\mathcal{L}}_{n}\left(P\middle|\widetilde{d}_{i,1}\right)=\left\{ \begin{array}{ll}\frac{1}{N-1}e^{-\frac{1}{2}\left(\frac{d_{i,1}-\sqrt{2}D_ {n}}{\sigma_{i,1}^{2}}\right)^{2}}&\text{if }\Delta d_{i,1}>\sqrt{2}D_{n}\\ \frac{1}{N-1}&\text{otherwise,}\end{array} \right. \tag{23}\]  
while the exact value of the likelihood function in center \(P\) is  
\[\mathcal{L}_{n}\left(P\middle|\widetilde{d}_{i,1}\right)=\frac{1}{N-1}e^{-\frac {1}{2}\left(\frac{d_{i,1}}{\sigma_{i,1}}\right)^{2}}. \tag{24}\]  
Thus, the HHT in a grid point \(P\) is the following:  
\[A(P)=\sum_{i=2}^{N}\mathcal{L}_{n}\left(P\middle|\widetilde{d}_{i,1}\right), \tag{25}\]and the upper bound of the HHT in this tile is  
\[\overline{A}(P)=\sum_{i=2}^{N}\overline{\mathcal{L}}_{n}\left(P\middle|\widetilde{ d}_{i,1}\right). \tag{26}\]  
Using the above results, Algorithm 1 describes the operation of the F-HHT. ```
Input:
1. measurements \(\widetilde{d}_{i,1},\sigma_{i,1},\ i=2,3,\ldots,N\)
2. sensor positions \((x_{i},y_{i},z_{i}),\ i=1,2,\ldots,N\)
3. desired grid size \(D_{final}\) Initialization:
4. create a grid over the search area \(\mathcal{S}\) with grid size \(D_{1}\)
5. add each tile to the PromisingList
6. \(GridSize:=D_{1}\)
7. set the minimal required value \(A_{maxmin}\) (e.g., to \(\frac{2}{N-1}\)) Iteration:
8. while \(GridSize>D_{final}/2\)
9. calculate \(A(P)\) and \(\overline{A}(P)\) for each tile in PromisingList
10. \(bestA:=\max(A(P))\)
11. \(bestP:=\underset{p}{\operatorname{argmax}}(A(P))\)
12. pruning 1: remove tiles with \(\overline{A}(P)<A_{maxmin}\)
13. pruning 2: remove tiles with \(\overline{A}(P)<bestA\)
14. \(GridSize:=GridSize/2\)
15. replace each tile with four smaller tiles of size \(GridSize\)
16. end Output: bestP
```  
**Algorithm 1** Operation of the F-HHT  
The inputs of F-HHT are the measurements, sensor positions, and the size of the grid on which the maximum is searched for (see steps 1-3). In step 4 of the initialization phase, the initial coarse grid is created. Then, in step 5, each tile is added to the list of promising tiles (PromisingList). The initial coarse grid size is set in step 6. The value \(A_{maxmin}\) in step 7 is the minimum required value of \(\max A(P)\), i.e., \(\max(A(P))>A_{maxmin}\). In a conservative solution \(A_{maxmin}=\frac{2}{N-1}\), since at least two hyperbolas intersect at the source position. When multiple sensors are available and there is a priori knowledge of the minimum number \(R\) of reliable measurements, then \(A_{maxmin}\) may be set close to, but below, \(\frac{R}{N-1}\). The iteration is performed between steps 8 and 16. In step 9, the F-HHT algorithm calculates the exact HHT (25) and the upper bound (26) for each promising tile center. The current highest exact value is \(bestA\) (step 10) at positions \(bestP\) (step 11). A tile with center \(P\) is non-promising (the maximum cannot be inside the tile) if \(\overline{A}(P)<A_{maxmin}\) or \(\overline{A}(P)<bestA\). In the pruning steps 12 and 13, the non-promising tiles are removed from the further search. In step 15, promising tiles are replaced by four new promising tiles with smaller size (see Figure 3). When the required grid size is reached, the center(s) of tile(s) with the highest HHT (\(bestP\)) form the output. Note that the output \(bestP\) may contain multiple grid points (tile centers). In such cases, the location estimate can be calculated as the mean of \(bestP\).