Joint Ultra-wideband Characterization of Azimuth, Elevation and Time of Arrival with Toric Arrays

In this paper, we present an analytical framework for the joint characterization of the 3D direction of arrival (DoA), i.e., azimuth and elevation components, and time of arrival (ToA) in multipath environments. The analytical framework is based on the use of nearly frequency-invariant beamformers (FIB) formed by toric arrays. The frequency response of the toric array is expanded as a series of phase modes, which leads to azimuth-time and elevation-time diagrams from which the 3D DoA and the ToA of the incoming waves can be extracted over a wide bandwidth. Firstly, we discuss some practical considerations, advantages and limitations of using the analytical method. Subsequently, we perform a parametric study to analyze the influence of the method parameters on the quality of the estimation. The method is tested in single-path and multipath mm-wave environments over a large bandwidth. The results show that the proposed method improves the quality of the estimation, i.e., decreases the level of the artifacts, compared to other state-of-art FIB approaches based on the use of single/concentric circular and elliptical arrays.


I. INTRODUCTION
I MPROVING the performance of wireless links requires a proper characterization and knowledge of multiple channel parameters: direction of arrival (DoA), time of arrival (ToA), delay spread, path loss, and K factor, among others [1], [2].With knowledge of the channel parameters, different scenarios can be effectively distinguished [3], even recreated and emulated through the use of post-processing techniques based on the creation/removal of reflections in the communication channel [4], [5].This is fundamental in today's interconnected society, considering the huge variety of propagation environments associated with different communication scenarios: This work has been supported by grant PID2020-112545RB-C54 funded by MCIN/AEI/10.13039/501100011033.It has also been supported by grants PDC2022-133900-I00, PDC2023-145862-I00, TED2021-129938B-I00 and TED2021-131699B-I00 funded by MCIN/AEI/10.13039/501100011033and by the European Union NextGenerationEU/PRTR; and by Universidad de Granada through grant PPJIB2022-05; and in part by the Predoctoral Grant FPU19/01251 and FPU19/04085.(Corresponding author: Alejandro Ramírez-Arroyo.) Alejandro Ramírez-Arroyo, Pablo Padilla and Juan F. Valenzuela-Valdés are with the Department of Signal Theory, Telematics and Communications, Research Centre for Information and Communication Technologies (CITIC-UGR), Universidad de Granada, 18071 Granada, Spain (e-mail: alera@ugr.es;pablopadilla@ugr.es;juanvalenzuela@ugr.es).
Among all the channel parameters, DoA and ToA are among the most noteworthy.A joint estimation of DoA and ToA is essential because current communication scenarios suffer from temporal variations of the physical properties of the channel.Accurate DoA and ToA information can account for the changes occurring in the channel.In the sub-6 GHz regime, DoA and ToA estimation techniques have been widely explored in the literature [17].Most of the applied techniques are of narrowband nature [18]- [24], which generally limits their use for modern broadband applications in the mm-wave frequency range.As a consequence, new efficient wideband alternatives are being sought that feature DoA and ToA simultaneously.For instance, narrowband approximations can be extended towards wideband applications through the use of time-delay beamformers [25].These approaches propose a frequency-dependent phase shift, which is applied on a narrowband decomposition of the wideband communication channel.Thus, the combination of narrowband nature approximations and frequency-dependent phase shifters leads to wideband approaches.One step beyond, the whole wideband channel can be considered through the design and implementation of nearly frequency-invariant beamformers (FIB), which has proven to be a suitable alternative for (ultra)wideband DoA estimation [26]- [28].The objective of FIB is to parameterize the array coefficients so that the spectral and spatial dependences can be treated independently [29].Previous implementations of FIB have ranged from the use of one-dimensional (1D) arrays [30] to two-dimensional (2D) configurations based on the use of circular arrays [31]- [34].Recently, the use of FIB was extended by the authors to include elliptical geometries [35], [36].This is a generalization of previous approaches based on circular geometries, as circular and linear arrays are subcases included in the more general elliptical arrays.However, many of the DoA and ToA estimation methods are only capable to accurately estimate the time of arrival and one of the two spatial DoA components, i.e., either azimuth or elevation, but not both at the same time [31]- [33], [35], [36].In this regard, it is not so common to find analytical approaches that accurately estimate azimuth, elevation and time of arrival in one go.
In this paper, we present a novel estimation method for joint 3D DoA (azimuth, elevation) and ToA characterization.The technique is based on the use of nearly frequency-invariant toric arrays.Following a similar rationale than in [31]- [38], the multipath frequency response, acquired in the different spatial points that are part of the toric array, can be expanded as a series of phase modes and a preselected frequency-dependent filter.This leads to diagrams in the azimuth-time (AoA-ToA) and elevation-time (EoA-ToA) domains, from which the direction and time of arrival of the incoming waves can be accurately estimated in a wideband range of frequencies in either single-path or multipath environments.This is a remarkable feature of the proposed method, as a joint 3D-DoA and ToA characterization is rarely found in the literature.Moreover, the accuracy in the joint DoA and ToA estimation is improved with respect to previous state-of-art frequencyinvariant beamformers.This is essentially due to the geometry of the torus.The geometrical disposal of spatial samples in a toric array efficiently exploits the arrangement for an optimal estimation.Within a torus, multiple concentric circular arrays can be defined in horizontal planes, thus improving the estimation of the azimuth angle.Then, the optimal vertical plane, which includes a single circular array, is selected for a precise estimation of the elevation angle and time of arrival.Finally, it is worth mentioning that the proposed technique works efficiently in under-sampling conditions, which can be used to reduce processing time and the number of samples employed.Thus, the main contributions of this work are summarized as: (i) Development of the theoretical framework for the use of frequency-invariant beamformers in three-dimensional geometries.The distributions of samples in threedimensional spaces provided by toric arrays allow the development of expressions that lead to a three-dimensional characterization of the direction of arrival (azimuth and elevation) and time of arrival.(ii) Method validation through simulation in the mm-wave frequency range.This approach can be employed in wideband single-path and multi-path channels by decoupling the spatial and temporal domains given the channel frequency response of each spatial sample.(iii) Parametric analysis of the toric geometry.A study of the geometrical parameters of the torus given the above framework is performed to define the optimal region of operation of the method.
These contributions open up the possibility of angular and temporal characterization of a 3D multipath scenario.This is achieved by means of the frequency-invariant beamforming framework proposed for toric arrangements throughout the paper.Moreover, as it will be detailed later, the proposed framework improves the quality of the DoA and ToA estimation compared to other state-of-the-art FIB approaches based on the use of circular and elliptical arrays.The document is organized as follows.Section II presents the mathematical framework for the joint characterization of the azimuth, elevation and time of arrival.It also discusses some practical considerations, advantages and limitations of the method.Section III illustrates some numerical examples, including single-path and multipath scenarios, in order to validate the proposed theoretical framework.A parametric study on how the main involved parameters affect the performance of the method is also carried out.Finally, general conclusions are drawn in Section IV.

II. THEORETICAL FRAMEWORK
A torus is defined as a closed surface formed by the Cartesian product of two circles.Parametrically, it can be defined as: where R is the distance from the torus center to the tube center, and ρ is the tube radius.φ is the azimuth angle, which represents the rotation around the axis of revolution.θ is the elevation angle, i.e., the rotation angle around the tube.R > ρ is considered in order to ensure a ring torus shape.
Due to the nature of this geometry, circles can be generated in the XY plane given a specific z height.Since z(φ, θ) only depends on θ for any value of φ, θ can be directly substituted in x(φ, θ) and y(φ, θ) as sin(θ) = sin(arccos(−z/ρ)) = 1 − (−z/ρ) 2 ).This expression can be further simplified, leading (1) to: (2) An example of circumference in this XY plane is illustrated in Fig. 1(a) when z = 0. Therefore, several circumferences can be arranged given the horizontal XY plane for several z values.Equivalently, two circumferences are found in the vertical plane formed by the Z-axis and a given angle φ.Fig. 1(b) shows this plane when φ = φ l , which represents a smart selection of the plane given the φ-axis.Mathematically, the Φ l Z plane is defined by the line that lies in the Cartesian z-axis with unit direction vector u = ẑ, and the line that crosses the origin {0,0,0} and the coordinate {cos(φ l ), sin(φ l ), 0} with unit direction vector v = cos(φ l )x+sin(φ l )ŷ.Thus, the Φ l Z plane is defined by the plane equation sin(φ l )x − cos(φ l )y = 0. Some examples of these arrangements are illustrated in Fig. 2.These distributions will be analyzed in detail in later sections.These properties of the torus to define circumferences in several planes will be used to accurately estimate the 3D DoA (azimuth φ and elevation θ), as well as the ToA τ in multipath environments.Now, let us assume an incident spherical wave l characterized by the azimuth angle φ l , elevation θ l and time of arrival τ l .This signal impinges on a torus formed by P samples on each circumference arranged in the XY plane, and P samples on each circumference of the Φ l Z plane, for a total of P 2 spatial samples.The estimation of the three previous parameters is performed in two steps: (i) The first one takes advantage of the circumferences located in the XY plane to make an accurate estimation of φ l .(ii) The second step makes use of the previous value of φ l to accurately estimate θ l and τ l given the Φ l Z plane.

A. Azimuth of Arrival (AoA)
In the first step, the frequency response at the center of the circumferences lying in the XY plane is expressed as where κ l stands for the complex amplitude for the l-th wave and τ l,z for the delay of the l-th wave for a given height z.The term τ l,z can be expressed as τ l,z = τ l + τ z , where τ l is the delay characterized at the center of the torus (z = 0) and τ z is an additional delay introduced by the height z.According to the torus geometry, the largest difference between τ z values is given by 2ρ/c, which is the diameter of the torus tube divided by the wave propagation speed, i.e., the speed of light c.This is the case for θ l = 0°and θ l = 180°.For θ l = 90°, the incident wave impinges on the different circles at the same time, so the difference between values of τ z becomes zero.Generally, it is satisfied that τ l ≫ 2ρ/c, so we can assume that τ l,z ≈ τ l .This effect generates a slight temporal spread in the time domain estimation, which will be analyzed in the following sections.Regardless of this fact, note that there is no negative effect on the estimation since the ToA is not yet estimated in this step.
For each circle at a height z, H p,l,z (f ) reads as the channel frequency response acquired at the p-th sample.Since the spatial samples are uniformly distributed around the position of H l,z (f ), H p,l,z (f ) is given by where d l,z = c τ l,z , and d p,l,z is the distance traveled by the wave to reach the p-th spatial sample.Therefore, the term is the attenuation factor given the distance between the torus center at a height z, and the p-th sample for a path loss exponent γ.The complex exponential term stands for the phase shift introduced by the additional distance to the center of the torus, where According to the torus geometry, the previous term can be expanded as where φ p,z is the azimuth angle for the p-th sample located at a height z, which is evenly spaced in the angular domain φ ∈ [0, 2π).The term r z stands for the radius of the circles deployed at different heights.From (2), it is defined as Through the expansion in Taylor series, ( 6) can be approximated as which allows us to simplify (4) to As a starting point, a first approach considers sin(θ l ) = 1 by fixing elevation incidence angle θ l = 90°.This approximation is fundamental since it will later allow the decoupling of the azimuth angular term φ l from the frequency term f through the Jacobi-Anger identity.Note that this approach is feasible because FIBs that have been developed to be robust over multiple elevation angles will be applied in later steps [37], while the accurate estimation of the elevation angle θ l is performed in Section II.B.Additionally, d l,z /d p,l,z ≈ 1 for d l,z ≫ ∆d p,l,z , i.e. the wave source is far away from the toric array (plane wave approximation).Note that this approach is necessary for the development of the theoretical framework throughout this section.However, numerical simulations and experiments performed with circular and elliptical arrays have shown that the characterization of the propagation channel is feasible, even in the near-field region, i.e. considering spherical wave propagation [35], [36].After these two considerations, (9) can be simplified to As observed in (10), phase and frequency components are linked in the complex exponential term.The Jacobi-Anger identity, given by [39] allows us to decouple phase (φ l ) and frequency (f ) terms, extending (10) to (12) The former step opens up the possibility of a solution based on nearly frequency-independent beamformers (FIBs).Note that J n,z (•) is the Bessel function of the first kind of order n for the circular array with radius r z .
The former expression can be projected onto basis functions of the form e jmφp,z that excite the frequency response H p,l,z (f ) when φ p,z = φ l .This method is known as phasemode expansion [31], [38], which provides the phase-mode domain H m,l,z (f ): e j(m−n)φp,z P .
(13) Although the latter expression may seem impractical from a computational perspective, note that the rightmost term (1/P ) P −1 p=0 e j(m−n)φp,z becomes zero when n = m.Otherwise, when n = m this term is one, leading to: Thus, this expression contains the phase-and frequencydecoupled terms, e jmφ l and e j2πf τ l,z , the latter implicitly included in H l,z (f ).Additionally, this excited frequency response H m,l,z (f ) includes the frequency-dependent component, which can be eliminated by the optimal choice of an inverse filter W m,z (f ).Mathematically, the phase-mode response H m,l,z (f ) can be calculated as with where J ′ m,z (•) is the first derivative of J m,z (•).This filter has demonstrated to be an optimal choice since J ′ m,z (•) avoid deep nulls in W m,z (f ), providing larger bandwidth for the estimation and performing an accurate estimation of the azimuth angle φ l for several angles θ l [37].
In order to improve the estimation, the phase-mode average domain between all the circles at different heights z is calculated, being the phase-mode expansion The distribution of the P circles in the XY plane is shown in Fig. 2.These circles are distributed in the range of heights Previous work has shown that the use of multiple arrays significantly improves the estimation of DoA and ToA [32], [35], [36].H l (f ) stands for the frequency response at the origin.Additionally, note that weights w i with i w i = 1 might be included in (17) to prioritize a specific subset of circles within the toric array in the estimation process.Although the use of weights is not mandatory in this development, as the toric geometry ensures a minimum radius r z , these weights might be necessary in other threedimensional geometries where one of the circle estimations degrades the overall characterization.
In multipath environments (l > 1), the phase-mode response for every wave impinging the toric array is the sum of the individual contribution from each wave, which leads to The 2-D Discrete Fourier Transform (DFT) of the former expression provides the joint azimuth DoA and ToA estimation as: where K stands for the number of frequency samples and frequency spacing f s = B/K with B being the signal bandwidth and M is the total number of considered phase modes.The term f min is the lowest frequency considered in the bandwidth and m and k are the indexes that characterize M phase modes and K frequency samples, respectively.The ToA τ l is extracted through the frequency f obtained in e j2πf τ l , while φ l is related to the phase-mode domain in the complex exponential e jmφ l [see (15)].Consequently, the time domain resolution and the maximum observable time are given by 1/B and (K − 1)/B.Similarly, the angular resolution is given by 2π/M in the phase domain.The previous expression can be efficiently calculated through the Fast Fourier Transform (FFT) algorithm.Note that, as a 2D-DFT based method, the resolution in the angular and time domains is determined by B and M .Therefore, it is advised to consider ultra-wideband channels and a high number of phase modes to maximize the resolution, and thereby mitigate the off-grid phenomena.

B. Elevation of Arrival (EoA) and Time of Arrival (ToA)
The second step takes advantage of the previous estimation of the azimuth φ l to accurately estimate the elevation angle θ l and the time of arrival τ l .Basically, the same P 2 spatial samples previously defined to generate P circles with P samples in the XY plane can be used to form circles in the planes defined by the φ−Z axes.These P circles, also with P samples per circle, are those contained in perpendicular cuts to the tube [see Fig. 1(b)].
Similarly to (3), the frequency response at the center of the tube, given an angle φ l , is where η l is the complex attenuation, and τ l,φ l is the delay of the wave l from the source to the tube center at φ l .Therefore, the frequency response at the p-th sample reads with After Taylor's series expansion (8), we arrive to Note that, by choosing the ring coincident with the φ l estimation, we ensure to align the wave incidence plane with the distribution of the P tube samples.Thus, the sine dependent variable term that appeared in eqs.( 6), ( 8) and ( 9), becomes constant and is removed from eq. ( 23).Visually, the chosen circular array laying in the Φ l Z plane is illustrated in Fig. 2.This fact ensures the correct estimation of θ l and τ l,φ l as the optimal estimation is provided when the plane of incidence coincides with the plane where the spatial sampling lies [35].
Therefore, considering a negligible attenuation between the edge and the center of the tube, H p,l,φ l (f ) is given by The former expression is similar to the one proposed in (10), with the difference of: (i) considering a radius ρ; (ii) the cosine term depends on the elevation θ l ; and (iii) the considered center is that of the tube and not that of the torus.Hence, a development in the form of basis functions e jmθ p,φ l , and the choice of an optimal W m,φ l (f ) filter, similar to that in (11)- (16), results in The summation of the multiple incident waves (18) leads to H m,φ l (f ).Finally, the application of the 2-D FFT to H m,φ l (f ) provides the joint angular-time domain for the DoA elevation angle θ and the ToA at the center of the tube τ φ l , i.e., H(θ, τ φ l ).According to the torus geometry, the delay at the center of the tube, τ l,φ l , is related to the delay at the center of the torus, τ l , as follows: This last step completely characterizes the DoA and ToA of the l-th wave by extracting the angles φ l , θ l , and the time τ l .

C. Considerations of the method
Subsections II.A and II.B have presented the theoretical framework for the 3D characterization of the communication channel.Although the method has been shown to work correctly in 2D scenarios with several geometries, due to some simplifications that are carried out, it is necessary to perform a correct assignment of the considered parameters.These approximations involve the appearance of artifacts, i.e. spectral contributions in H(φ, τ ) and H(θ, τ φ l ) that may mislead with the real incident wave.As long as these artifacts are controlled, the characterization can be properly performed.For this purpose, we define a metric ∆, which indicates the difference between the power of the incident signal and the largest artifacts.Therefore, ∆ shows the array response of the toric geometry when a wave with certain parameters (φ l , θ l , τ l ) impinges on the array given the joint angular-time domains H(φ, τ ) and H(θ, τ φ l ).If the value of ∆ is greater than 0 dB, the incident signal is distinguishable from the artifacts [35, Fig. 3].Thus, large ∆ values maximize the dynamic range when characterizing the communication channel.
Concerning the number of spatial samples P , it must be chosen in such a way that we ensure compliance with the Nyquist spatial sampling theorem, i.e., a separation between samples less than half wavelength.Otherwise, spatial aliasing appears, which translates into an increase in the level of the artifacts.Fundamentally, this depends on the radii considered for the circles in the different planes of the torus.In Section II.A, the boundary is given by the outer circle when z = 0, where r z = R+ρ, while in Section II.B, the considered radius is directly ρ.
Given the maximum value of r z , we can ensure that the sampling theorem is satisfied on the torus with P equispaced samples per circle if Regarding the number of phase modes M , it can be demonstrated that J n (•) ≈ 0 for a sufficiently large order n [33], [39].Thus, M must be chosen so that it is lower than the above threshold.Otherwise, the denominator in ( 16) would tend to zero, generating numerical instabilities in the method.Logically, the greater the number of filters, the greater the required computational time.The number of raw W m,z (f ) filters needed is M × P in the first step.However, by taking advantage of the reflection (mirror) symmetry of the torus on the z-axis, the number of filters can be decreased to M × (P/2 + 1).Also, by taking advantage of the symmetry of the Bessel functions for positive and negative indexes m [J −m (x) = (−1) m J m (x)], the total number of filters is reduced to (M/2 + 1) × (P/2 + 1).In the second step, only a single circle is considered, thus, M/2 + 1 filters are needed in this case.In total, (M/2 + 1) × (P/2 + 2) filters are required for a joint characterization of the AoA, EoA and ToA.

III. VALIDATION OF THE METHOD
This Section discusses and validates several aspects of the method presented in Section II.Specifically, a parametric study shows the main implications of the different parameters involved in the joint characterization of the AoA, EoA and ToA.According to the torus geometry, several aspect ratios, i.e., R/ρ, and sizes are considered to show the good performance of the method for diverse cases.Later, the joint estimation is validated through simulations at several frequency ranges for different scenarios based on both single-path and multipath scenarios.

A. Parametric analysis
Regarding the physical geometry of the torus, R, ρ and P are the main parameters to be considered.As previously stated, these are directly related to the Nyquist spatial sampling theorem [see (27)].In order to analyze the effect of fulfilling this theorem, it is considered a case with a single-path scenario impinging the toric array with τ l = 15 ns (d l = 4.5 m), and φ l = 180°for several values of θ l .The torus size is given by R = 0.5 m and ρ = 0.25 m, with a frequency range that goes from 28 GHz to 32 GHz (B = 4 GHz) and K = 200 frequency samples.M = 300 phase modes are considered.Given R, ρ, and f max = 32 GHz (λ min = 9.37 mm), the Nyquist theorem is satisfied even for outer circles (r z = 0.75 m) if P > 1006.For 336 < P < 1006, it is exclusively satisfied for some inner circles.If P < 336, it is not fulfilled even for the inner circles (r z = 0.25 m).
For the given parameters, Fig. 3 illustrates the value of the metric ∆ when the number of samples P is varied.Note that a higher value of ∆ implies a better DoA and ToA estimation, as the level of the artifacts is reduced.For the case θ l = 90°, the quality of the estimation is optimal due to the approximation performed in (10).However, even non-coincident waves in the elevation plane (θ l = 90°) provide good results for ∆.Concerning the number of spatial samples P , a flat behavior is observed in Fig. 3 when the Nyquist theorem is fulfilled (P > 1006).This implies that the estimation will not improve even if the number of samples is increased.In an intermediate region (336 < P ≤ 1006), even though sampling is not performed correctly for the outer circles, a flat behavior is still observed until the sampling theorem is no longer satisfied for about half of the circles, i.e., r z = R.This fact contrasts with estimation on single arrays, where not complying with the sampling theorem dramatically increases the artifacts [35].The geometry of the torus, which supports estimation from concentric circles [see (17)], allows the number of P sensors to be reduced below the sampling theorem, compensating for this decrease with the larger number of arrays deployed.When P is further reduced (P < 336), the value of ∆ begins to significantly decay.In any case, note that despite the presence of artifacts, ∆ > 0 dB, which makes the incident wave distinguishable from the artifacts even for a reduced number of samples.
In a second study, the influence of the number of phase modes M in the DoA and ToA estimation is analyzed.All the parameters remain the same as in the previous experiment, except from the torus geometry.The parameters of the torus are: R = 0.75 m, ρ = 0.25 m and P = 1440.Fig. 4 shows the value of ∆ when varying the number of considered phase modes M .Above a certain threshold, approximately M > 700 in this case, J n (•) ≈ 0. This causes numerical instabilities in the computation of the filter W m,z (f ) [eq. ( 16)] as the denominator approaches zero.As a result, the quality of the estimation significantly degrades.In the range 300 < M < 600, an almost flat region is found, showing that, regardless of the number M chosen, the estimation performs correctly.Finally, if M is even more decreased, the estimation tends to degrade, besides decreasing the angular resolution (see Sec. II.A).Ideally, one should operate in the nearly flat region where ∆ approaches the maximum.This region depends exclusively on the 2πf r z /c argument of the Bessel function, so a prior analysis of the optimal working regions is essential.
From a single AoA/EoA/ToA domain perspective, the array gain of the toric geometry can be obtained by averaging a single dimension of the phase-mode H m,l (f ) [see (17)] or H m,l,φ l (f ) [see (25)], given the pairs {φ l , τ l } or {θ l , τ l,φ l }, respectively.Fig. 5 shows the normalized array gain for several φ l , θ l and τ l values given the AoA/EoA/ToA domains.The high directivity, whether in angular or temporal domain, is illustrated across the three domains, which will enable accurate simultaneous characterization of azimuth, elevation and time of arrival.Finally, to demonstrate the potential of the toric geometry, the metric ∆ is compared with previous two-dimensional geometrical approaches that make use of frequency-independent beamformers and phase-mode transformations.For this purpose, we consider the same incident wave (φ l and τ l ) as in the previous experiments, and the same frequency range, B and K. Fig. 6 shows the value of ∆, given the azimuth (AoA) -time (ToA) domain, for different elevation angles θ l given four different geometries: (i) elliptical array, (ii) rotated elliptical array, (iii) circular array and (iv) toric array.Therefore, Fig. 6 compares the performance of FIBs combined with phase-mode transformations for several arrangements.Specifically, the theoretical framework for elliptical arrangements has been developed in our previous works [35], [36], while the framework for circular arrangements is based on [32], [37].The number of phase modes is fixed to M = 300.Both elliptical arrays have a semi-major axis of 0.5 m, an eccentricity of 0.7, and P = 720.Additionally, the rotated ellipse includes an angle rotation of 45°.The single circular array has a radius of 0.5 m, and P = 720.The toric array parameters are R = 0.5 m, ρ = 0.25 m and P = 720.It can be appreciated that the toric array estimation outperforms all other geometries.This is due to the joint estimation based on the use of the concentric circular arrays that are discerned in the toroid geometry itself.Note that this sample distribution will later allow us to perform the accurate estimation of the elevation angle θ l , while in the other geometries, it is not possible due to the two-dimensional distribution of the samples.The application of the two-stage method might imply the presence of propagation errors if θ l is calculated based on an inaccurate characterization of φ l , or vice versa when φ l is initially calculated with an unknown θ l .Nevertheless, Fig. 6 illustrates the robustness of the method given the characterization for several angles.
In summary, the proposed method is shown to improve the quality of the estimation by decreasing the level of the artifacts compared to other geometries, such as single circles or ellipses, even enabling the estimation with under-sampling conditions due to the use of several arrays, i.e., several FIBs, simultaneously.

B. Single-path Characterization
Once some of the key parameters for the method have been analyzed, this Section illustrates some examples of the 3D joint characterization of the AoA (φ l ), EoA (θ l ), and ToA (τ l ).Let us assume a frequency range f ∈ [58, 62] GHz with B = 4 GHz and K = 200 frequency samples.According to Section II.A, for R = 0.25 m and ρ = 0.125 m, P = 720 spatial samples and M = 300 phase modes are enough to reduce The level of the artifacts is below -30 dB, so they do not appear in the figure.This effect is obtained due to the joint estimation of the multiple concentric circles contained in the torus.
Given the azimuth estimation, the circumference lying on this angle φ l is chosen.Fig. 7(b) illustrates the elevation-time domain.Note that the value of ∆ is maximum at the EoA (θ l = 90°) and ToA (τ l,φ l = 19.25 ns).According to eq. ( 26), the exact time of arrival should be 19.17ns, appearing this difference due to the temporal resolution 1/B.Nonetheless, this value is within the resolution range of the method, being a correct τ l,φ l estimation.In Fig. 7(b), the level of the artifacts is below -23 dB.
In a second experiment for the single-path scenario, we consider a wave with φ l = 160°, θ l = 110°, and τ l = 30 ns.The additional parameters remain identical to the previous experiment.Note that the temporal dispersion explained in Section II.A due to an elevation θ l = 90°is not appreciable because of the high working frequencies (f c = 60 GHz), since the physical size of the array decreases notably, being τ l,z ≈ τ l .Therefore, the estimation τ l performed in the first step can be considered as reliable taking into account high frequencies in the mm-wave band with small array sizes where τ l ≫ 2ρ/c.The two previous experiments have shown the ability to correctly identify waves in the three-dimensional space and in the mm-wave frequency range.

C. Multipath Characterization
The previous section has validated the performance of the FIB for the 3D channel characterization in a single-path scenario.According to (18), for a multipath scenario (l > 1), the spectral response of the beamformer is directly the sum of the responses of each incident wave.However, in practice, the amplitude of the spectral response is actually dependent on the elevation θ l due to the assumption made in (10).This fact was corroborated in Fig. 6 with differences of up to 30 dB when going from θ l = 90°to θ l = 130°.This implies that, under the assumption of a set of L incident waves with similar θ l , the DoA-ToA domain provides a correct characterization of the scenario.However, if this set of waves has θ l different from each other, the waves with θ l close to 90°may fade the rest of the incident waves.Therefore, estimation in multipath environments with several θ l values, i.e., the most general case, can be performed based on a subtraction strategy of the previ-This work has been accepted for publication in IEEE Transactions on Wireless Communications.DOI: 10.1109/TWC.2024.3377539ously estimated waves.Given the estimation of the l-th path, the influence of the previous l − 1 paths can be removed as L n=l H p,l,z (f ) = L n=1 H p,l,z (f ) − l−1 n=1 H p,l,z (f ).If the estimation is performed correctly, the effect of the substracted l − 1 waves is removed from H(φ, τ ) and H(θ, τ φ l ).In case that the subtraction is not effective due to the temporal (1/B seconds) and angular (2π/M radians) resolution values, one can proceed to reconstruct H p,z (f ) as a discretized set of waves based on a low-complexity search with expectationmaximization algorithms such as Space Alternating Generalized Expectation-maximization (SAGE) [34], [40].Thus, by employing this approach, the off-grid phenomena can be mitigated.
To summarize the rationale previously described, Algorithm 1 presents a pseudocode that illustrates the main steps to be followed in a general multipath environment.The method, originally based on DoA and ToA estimation for circular arrays, has been generalized to toric arrays where the geometry of the torus has been exploited to define circles in multiple planes.Through an appropriate choice of the parameters involved in the method, good performance of the characterization will be obtained.
In order to validate the multipath characterization, a scenario with L = 3 waves is simulated.

Fig. 1 .
Fig. 1.Torus geometry: Representation of the circles arranged in (a) the XY plane, and (b) Φ l Z plane.

Fig. 2 .
Fig. 2. Planes where the circular arrays are placed on the torus.P circular arrays can be found in the XY plane given several heights z (left panel).For visualization purposes, only five circles are marked in blue in the XY plane.A single circular array is located in the Φ l Z plane (right panel).The P circular arrays in the XY plane and the single array in the Φ l Z plane allow us to estimate the AoA and the pair {EoA, ToA}, respectively.

Fig. 3 .
Fig. 3. Relation between the power of the estimated wave and artifacts (∆) as a function of the number of spatial samples P for several incident elevation angles θ l .Parameters of the considered scenario: R = 0.5 m, ρ = 0.25 m, φ l = 180°and τ l = 15 ns.Non-compliance with the spatial sampling theorem increases the artifacts, thus decreasing the value of ∆ in the estimation.

Fig. 4 .
Fig. 4. Relation between the power of the estimated wave and artifacts (∆) as a function of the number phase modes M for several elevation angles θ l .Parameters of the considered scenario: R = 0.75 m, ρ = 0.25 m, φ l = 180°a nd τ l = 15 ns.
Figs. 8(a) and 8(b) illustrate the azimuth-time and elevation-time domains, respectively, for the joint estimation of this incident wave.As observed, the maximum value of ∆ coincides with the AoA, EoA and ToA positions in both figures.Artifacts amplitude appears 20.3 dB and 18.4 dB below the correct estimation in Figs.8(a) and 8(b), respectively.