Spherical Fourier-Transform-Based Real-TimeNear-Field Shaping and Focusing in Beyond-5G Networks

For ultra-reliable high-data-rate communication, the beyond fifth generation (B5G) and the sixth generation (6G) wireless networks will heavily rely on beamforming, with mobile users often located in the radiative near-field of large antenna systems. Therefore, a novel approach to shape both the amplitude and phase of the electric near-field of any general antenna array topology is presented. Leveraging on the active element patterns generated by each antenna port, the beam synthesis capabilities of the array are exploited through Fourier analysis and spherical mode expansions. As a proof-of-concept, two different arrays are synthesized from the same active antenna element. These arrays are used to obtain 2D near-field patterns with sharp edges and a 30 dB difference between the fields’ magnitudes inside and outside the target regions. Various validation and application examples demonstrate the full control of the radiation in every direction, yielding optimal performance for the users in the focal zones, while significantly improving the management of the power density outside of them. Moreover, the advocated algorithm is very efficient, allowing for a fast, real-time modification and shaping of the array’s radiative near-field.


Introduction
Key performance indicators for the evolved fifth generation (E5G) and upcoming sixth generation (6G) mobile wireless communications networks promote continuous connection availability, strong reliability, huge device density and low air interface latency [1]. To unleash the full potential of E5G and 6G for indoor applications, a holistic multi-disciplinary approach is required based on disruptive communication technologies and innovative beamforming architectures [2][3][4]. In particular, to support extremely high data rates and to improve link reliability, advanced beamforming networks operating in millimeter-Wave (mmWave) and TeraHertz (THz) frequency bands play a key role. In such cases, however, the near-field (Fresnel) distance can amount to several dozens of meters. Therefore, the proper assumptions should be adopted to analyze the system's performance, instead of the conventional far-field approach [5][6][7][8].
To enable antenna arrays to set up wireless communication in the Fresnel region, an appropriate focusing technique is essential to obtain properly shaped focal spots exhibiting high energy concentrations in the near-field region [9]. Moreover, to be able to modify the shape and position of focusing in real-time, the array's antenna elements' excitations require fast updates in terms of new amplitudes and phases. Therefore, the design of a practically feasible near-field focused array (NFFA) must fulfill stringent and sometimes conflicting requirements, often not occurring in standard base stations. On the one hand, the use of complex topologies or expensive materials must be prevented to reduce cost and footprint [10,11]. On the other hand, to achieve an improved wireless efficiency, This paper is organized as follows. Section 2 outlines the problem statement, details the far-to-near-field transformation and the subsequent shaping of the fields radiated by the array aperture's current distribution, based solely on the compressed active far-field data. Section 3 validates and demonstrates the technique through both numerical and full-wave experiments. Conclusions and future plans are summarized in Section 4.

Problem Statement
Consider the application scenario of Figure 1, where a number of users are distributed in the radiative near-field of an antenna array. The goal of this work is to determine the complex excitation currents required to shape the array's electric near-field E(r) in order to ensure user-specific and quasi interference-free operation.
Antenna element An ten na arr ay The array is conceptually represented in Figure 2a. Attached to its phase center is the origin of the global coordinate system O(x, y, z), which also coincides with the center of the array's circumscribing sphere of radius R. In the following sections, E(r) is shaped so that it corresponds to a prescribed target-field distribution T(R Tr ) over the surface of a sphere of radius R T ≥ R (Figure 2a), wherer denotes the radial unit vector. Therefore, we propose a novel algorithm that efficiently computes the required excitations of the array's antenna elements, using only their active radiation patterns as input.

Far-to-Near-Field Transformation
Consider the nth antenna element with its phase center indicated by the vector p n and attached to the origin of a local coordinate system O n (x n , y n , z n ) (Figure 2b). To obtain the nth active far-field radiation pattern, this element is excited via its port by a unit current I n = 1 A injected by an appropriate Norton equivalent source I g,n . All other antenna elements i = n are terminated by a 50 Ω impedance. This is schematically shown in Figure 3. Due to mutual coupling, the excitation gives rise to a current distribution j n (r n ), flowing over the whole volume V, i.e., on all antenna elements in the array, where r n is defined in the local coordinate system O n .
(a) Conceptual sketch of an array of volume V circumscribed by a sphere of radius R and with its phase center coinciding with the origin of the coordinate system O(x, y, z). The goal is to shape the array's emitted field so that it corresponds to a desired pattern on a sphere with radius R T . (b) Each nth array element, fed through the nth port, is attached to a local phase center, indicated by a vector p n , and coinciding with a local coordinate system's origin O n (x n , y n , z n ). The circumscribing sphere with regard to the local coordinate system has a radius R n . A current I n is injected into the nth port (see also Figure 3), generating a current distribution j n on the entire array.
arch 10, 2023 submitted to Sensors 5 of 16 I g,n Z n 50 Ω I n = 1 A Z i 50 Ω Figure 3. Excitation of the n th active far-field radiation pattern. The n th port is excited by a Norton equivalent source, formed by a current source I g,n in parallel with a 50 Ω resistance. The source is chosen such that a current I n = 1 A is injected into the n th antenna, represented by its input impedance Z n . All other ports for i ̸ = n are instead terminated by a 50 Ω resistor.
where the coefficients f m l,n are the result of a spherical Fourier transform of the n th currentnormalized and active far-field radiation pattern F n , as [39, eq (8)] As explained in [39], when a far-field is expressed in the spherical basis, as F = F θθ + F ϕφ , 129 both F θ and F ϕ are not continuous functions on the unit sphere Ω. In particular, they exhibit 130 discontinuities at the north and south poles. Therefore, to compute (6), we express the 131 radiation pattern in Cartesian components as F = F xx + F yŷ + F zẑ [40]. These radiation 132 patterns are, for instance, obtained via simulations or measurements.

133
Further, it is known that each F n is (quasi-)band-limited, meaning that it can be efficiently described using a L n -order spherical Fourier transform, as [39, eq. (9)] where ϵ(L n ) is the maximum absolute error in dBV stemming from (8). To minimize 134 ϵ(L n ), on the one hand, the f m l,n coefficients (6) are computed through a Q-order Lebedev 135 Figure 3. Excitation of the nth active far-field radiation pattern. The nth port is excited by a Norton equivalent source formed by a current source I g,n in parallel with a 50 Ω resistance. The source is chosen such that a current I n = 1 A is injected into the nth antenna, represented by its input impedance Z n . All other ports for i = n are instead terminated by a 50 Ω resistor.
The electric field E n (r n ), generated by the current density j n (r n ) and defined with respect to O n is then obtained by the electric field integral equation (EFIE) [35], which reads where ω is the angular frequency, and µ is the free-space permeability. In (1), we make use of a multipole expansion for the dyadic Green's function G [36] detailed in Appendix A and valid for r n = r n ≥ R n , i.e., where R n (Figure 2b) is the radius of the sphere circumscribing the current distribution j n (r n ) and thus the antenna array, measured with respect to O n . Further, k is the wavenumber, I represents the unit dyadic, andȲ m l is the complex conjugate of the l th order, m th degree scalar spherical harmonic [37]. The multipole function H m l in (2) is defined as where h (2) l is the l th order spherical Hankel function of the second kind ( [38], p. 437). Furthermore, throughout this Section, we make use of the following notation, Substituting (2) into (1) yields a compact multipole expansion of the electric field E n , valid for r n ≥ R n and thus including the array's radiative near-field, as where the coefficients f m l,n are the result of a spherical Fourier transform of the nth currentnormalized and active far-field radiation pattern F n , as ( [39], Equation (8)) F n (k) = − ωµ 4π I −kk ·ˆV e k·r n j n (r n ) dr n .
As explained in [39], when a far-field is expressed in the spherical basis as F = F θθ + F φφ , both F θ and F φ are not continuous functions on the unit sphere Ω. In particular, they exhibit discontinuities at the north and south poles. Therefore, to compute (6), we express the radiation pattern in Cartesian components as F = F xx + F yŷ + F zẑ [40]. These radiation patterns are, for instance, obtained via simulations or measurements.
Further, it is known that each F n is (quasi-)band-limited, meaning that it can be efficiently described using an L n -order spherical Fourier transform, as ( [39], Equation (9)) where (L n ) is the maximum absolute error in dBV stemming from (8). To minimize (L n ), on the one hand, the f m l,n coefficients (6) are computed through a Q-order Lebedev quadrature [41]. The value of Q yielding the best performance is obtained by rounding up L n to the closest available quadrature order. On the other hand, for a radiating structure with a radius R n , modes of order l > kR n provide little additional information about its radiation pattern and electric field [42]. As a rule of thumb, both the summations in (5) and (8) may therefore be truncated at L n = kR n .
Whereas a higher value of L n might reduce , a trade-off between accuracy and performance must be carefully considered. Empirically, it was verified that a value of maximally < 10 dBV is required to accurately model the device under test. To speed up this computation, the values that the functionȲ m l (r) in (6) assumes for all possible Lebedev quadrature nodes are stored during the set-up phase of the algorithm, amounting to circa 100 MB of data, and used for future simulations.

Phase Center Translation
Owing to linearity, the total electric field E(r) produced by the array is given by the weighted sum of the electric fields produced by all N antenna elements, as in where I n is the nth feeding current, and E n (r) is the electric field produced by exciting the nth antenna element, expressed with respect to the global phase center O. To perform a translation from the nth local to the global coordinate system, we make use of some simple vector algebra: where p n (Figure 2b) is the vector indicating the position of the nth local phase center in the global coordinate system. An addition theorem for (12) is easily derived from the classic theorem for spherical Hankel functions ( [43], Equation (6)). For where the translation operator T µm λl is defined as in which the symbol G indicates a Gaunt coefficient ( [44], Equation (A-2)) and where with j q being the qth order spherical Bessel function of the first kind. Plugging (13) into (12), we express the nth global electric field, valid for r ≥ R, as where (17) is a linear mapping that transforms the original set of coefficients f m l,n , determined up to order L n and with respect to the nth local coordinate system O n (x n , y n , z n ), into the translated coefficients f µ λ,n required up to order Λ n and with respect to the global coordinate system O(x, y, z).
To determine the optimal value of Λ n , we make use of the nth cumulative power spectrum Γ n , defined as ( [39], Equation (23)) which is a measure for the power radiated when exciting the nth antenna element. We exploit the property that, under translation of the antenna's phase center, the cumulative power spectrum remains constant. Theoretically, Λ n can be obtained from the following equivalence Numerically, we compute Λ n iteratively, since the new cumulative power spectrum Γ n (Λ n ) approaches Γ n (L n ) asymptotically for increasing Λ n , as shown in Figure 4. Simulations have shown that defining Λ n as the value for which Γ n (Λ n ) is at least 99% of Γ n (L n ) yields the best balance between computational time and accuracy.
. Typical sigmoid behavior of the cumulative power spectra ratio Γ n (Λ n )/Γ n (L n ) for increasing Λ n . A sufficient value for Λ n is reached when this ratio is at least 99%.

Near-Field Intensity Shaping
Assume that the antenna array is designed to emit an electric field with a specific polarizationû. We split each nth electric field E n as so that by combining (11), (16) and (21), the desired co-polar component of the array's electric field values, at a distance R T from the array's phase center, i.e., for r = R Tr , is given by where f µ,co λ,n = f µ λ,n ·û, and Λ = max (Λ n ). To achieve near-field shaping, we set the left-hand side of (23) equal to a target field T(R Tr ) = T(R Tr ) ·û. Subsequently, the currents I n that best approximate the desired pattern are obtained as It is important to note that the only restriction imposed by the addition theorem is that R T ≥ R, meaning that it is possible to obtain the values T(R Tr ) on the surface of a sphere located entirely in the array's radiative near-field. Combining (3) and (24) and performing a Λ-order multipoles expansion (8) For very large values of R T , the advocated method smoothly transitions into a more refined version of the traditional far-field beamforming technique. This is easily proven by calculating the denominator in (25) which is equivalent to a simple change of overall amplitude and phase reference of T(R Tr ).
In this way, for very large R T values, the advocated method translates to a far-field shaping technique. Finally, we collect the sought-after excitation currents I n into the unknown vector i and the target Fourier coefficients t µ λ into the known vector t. This way, we cast (25) into the following linear system of equations, where the elements of the system matrix M are the antenna element's coefficients f µ,co λ,n . As a result, M has dimensions M × N, where M = (Λ + 1) 2 is the number of harmonics involved in the computation, and N is the amount of antennas in the array. Before solving (27), however, we have to verify that the problem is well-conditioned. This is equivalent to the columns (or rows) of M being linearly independent and, as a consequence, a limited condition number κ(M).
The nth column of M collects the translated spherical Fourier spectrum's coefficients of the nth antenna element's active radiation pattern (17). Therefore, the elements of M are only a function of the original active radiation patterns and of their respective element positions p n . In general, for antenna arrays consisting of identical elements with sufficient inter-element spacing (i.e., several λ), the difference between all the active radiation patterns is practically negligible. In this scenario, κ(M) becomes solely a function of the array layout. In other words, by ensuring that the antennas are spaced sufficiently far from each other, the columns of M will be sufficiently linearly independent. For denser antenna arrays, it was empirically verified that a spacing of at least λ/3 yields values of κ(M) in the order of 10 3 ∼10 4 , while a spacing of λ/2 reduces it further to order 10 2 ∼10 3 . For very dense antenna arrays with a spacing smaller than λ/3, the problem of an increasing condition number is partially mitigated by the higher coupling between the antenna elements. Owing to the asymmetries in the array layout, a larger variation in the data is introduced, reducing the correlation between the rows of M. While in this case the set-up time would inevitably increase, the novel method is still applicable to very dense arrays owing to this set-up operation (i.e., obtaining the active far-field data and computing M) being performed only once per layout. However, special care should be taken when interpreting this superdirective solution, as it might be very sensitive to small changes in array geometry, excitations and operating frequency. On the other hand, it is known from the literature [45] that the value of Λ increases linearly with the array size. Since the value of M, and therefore the number of rows of M, scale quadratically with Λ, it is generally the case that M N, meaning that (27) corresponds to an overdetermined linear system. Since we can ensure that M has full column rank, we opt for a least-squares solution as where † indicates the conjugate transpose operation. Note that the construction of the system matrix M and the subsequent computation of its pseudo-inverse M + are part of the set-up process of the algorithm. Once M + has been determined for a specific array, the process of adjusting the target in both distance and shape and obtaining the necessary feeding currents simply boils down to two simple steps that can be performed in real time. First, the target vector t is determined via (25) by a Λ-order Lebedev numerical quadrature, requiring circa Λ samples. Its evaluation requires 2Λ multiplications and Λ additions. For all orders λ and degrees µ, the construction of t needs 3Λ(Λ + 1) 2 ≈ 3Λ 3 operations. Second, one matrix-vector multiplication is performed to evaluate (28). Since the pseudo-inverse matrix M + has dimensions N × M, while the vector t has dimensions M × 1, a total of 2N M operations are required. Since M = (Λ + 1) 2 , the total number of operations scales as 2NΛ 2 . Therefore, the near-field shaping procedure is performed in less than a second, rendering the proposed method suitable for real-time field-shaping applications.

Numerical Validation
To validate the novel method, we make use of an x-polarized half-wave dipole antenna of length 136 mm, 0.9 mm wire diameter and operating at the frequency of 1 GHz (λ = 30 cm) as the elementary antenna element in our array configurations. A planar array consisting of N = 197 such elements is then synthesized. The antennas are distributed inside a circular region of radius R = 4λ = 1.2 m, forming a grid with inter-element spacing of λ/2 = 15 cm. In this experiment, we choose the radius of the target sphere as R T = 6λ = 1.8 m. As this size amounts to approximately 75% of the array's aperture dimension, the shaping process takes place on a sphere in the radiative near-field of the structure.
The theory described in Section 2 dictates that each of the 197 active radiation patterns F n is to be determined, which represents quite a cumbersome task. However, owing to the rather large inter-element distance, the mutual coupling between neighboring antenna elements remains quite limited. Therefore, we may approximate all F n by the active radiation pattern F 0 of the element closest to the array's center, i.e., the one located at x = y = 0. This is obtained via a full-wave simulation performed by the frequency domain solver of CST Microwave Studio [46] and normalized according to Figure 3. These data are used as input for a Python [47] script implementing the advocated algorithm. The spherical Fourier transform of F 0 is computed via (6) with an accuracy of −10 dB, requiring a value of L 0 = 18. Using (17), the resulting f m l,0 coefficients are then translated to the 197 different antenna positions. According to (19), a value of Λ = max (Λ n ) = 30 is needed, corresponding to M = (Λ + 1) 2 = 961 harmonics involved in the numerical computation. This results in a 961 × 197 system matrix M, with a condition number κ(M) ≈ 3547. Subsequently, the target's spherical Fourier transform and its corresponding target coefficients t are determined via (25). Finally, the unknown vector i is obtained via (28). This process required only 230 ms on an Intel Core i7-8650U processor running at 1.90 GHz with 16 GB of memory. Two examples are shown in Figures 5 and 6, where the plots' labels are given by where θ ∈ 0, π 2 and φ ∈ [0, 2π]. In Figures 5a and 6a, the required excitation currents for all 197 antennas computed by the novel method for the two target fields are shown.
The current amplitudes produced by the advocated algorithm naturally exhibit a Gaussian-like distribution, with elements closer to the array's center receiving higher excitations. This is typical of tapering techniques [48], which are usually required in traditional array focusing methods as an additional step to reduce side-lobe levels. On the other hand, the current phases exhibit a checkerboard-like pattern, with neighboring elements having close to perfectly opposing phases. This is required to achieve radiationless interference between adjacent elements, allowing the remaining evanescent fields to carry shape features to the target surface [34]. The corresponding co-polar electric near-fields' values obtained for the upper hemisphere of the target surface are shown in Figures 5c and 6c. A very close resemblance to the target fields of Figures 5b and 6b, which are shapes with sharp edges, is obtained. While at first sight the phase profile is not overall flat, the novel method ensures the correct values over the areas of interest, i.e., where the shapes are defined. Note that the phase values are irrelevant in areas where the magnitudes of the field values are very low. Since the advocated algorithm operates in the Fourier domain, where a large amount of far-field data is replaced by a compact angular spectrum description, it requires only about 100 ms to import any new target shape and then determine the currents necessary to obtain it. As previously mentioned, our novel method is therefore applicable to real-time field-shaping.

Full-Wave Validation and Application
To further validate the method, the results of the synthesis using the advocated technique are compared to accurate full-wave simulations using CST Microwave Studio. To limit the simulation time in CST, we design a 7 × 7 planar array using the same dipole antenna as antenna element. The inter-element spacing is 0.6λ = 18 cm, for a total size of circa 121 × 108 cm 2 . The radius of the target sphere is now R T = 4λ = 1.2 m, corresponding to one array aperture, to guarantee that the shaping process happens on a sphere located in the radiative near-field. The target fields are shown in Figure 7b,f. As in the previous section, the f m l,0 coefficients pertaining to the center antenna element are processed by (17), yielding a value of Λ = max (Λ n ) = 22. This corresponds to M = (Λ + 1) 2 = 529 harmonics involved in the numerical computation, resulting in a 529 × 49 system matrix M, with a condition number κ(M) ≈ 3. Computing the required currents for this array took circa 160 ms. Figure 7a,e show the feed currents that are computed by means of the novel method in order to obtain the target fields. The resulting field values are shown in Figure 7c,g. We now also use the computed currents to feed the array in CST Microwave Studio and perform a full-wave simulation of the array's resulting fields. These are shown in Figure 7d,h. Finally, we return to the application scenario of Figure 1, where U users are located close to one another. To ensure a good coverage for each individual user, the currents (and resulting electric near-field) obtained from the experiment of Figure 7a are combined in order to focus the power separately at U different spots. Furthermore, through a proper implementation of the feeding architecture of the array system, this permits us to deliver U different data streams to the users. The outcome of such a scenario, with U = 4, is shown in Figure 8.
When comparing target patterns, numerical solutions and full-wave results, an excellent agreement is observed, in both amplitude and phase. Any small difference is due to the approximation made by employing the center element active radiation pattern to represent all the remaining elements. This choice is, however, justified by both the enormous amount of time saved in performing the simulations of the antenna patterns and by the very good match between the numerical results and the full-wave simulations. Note that compared to the previous section, there is a loss of resolution as the sharp edges of the target fields are less accurately reproduced. Nevertheless, this is merely due to the use of a smaller array aperture rather than to the algorithm's loss of accuracy.

Discussion on Computational Efficiency
All the methods referenced in this work, including the presented technique, require a one-time set-up operation. While an accurate description of this set up is not always provided by the respective authors, it can be estimated to be approximately the same across all algorithms. Therefore, we focus the following comparison on the time required for the actual electric near-field shaping process across different methods: • Smart skin holography [13]: While a quite good resolution is obtained in near-field shaping, the process required circa 20 min per shape. Again, this algorithm is also not suited to real-time applications. • Method of Moments (MoM) [30]: A radial slot array is synthesised with an in-house MoM code and tailored for a specific target shape. The amount of time usually required to fill an MoM matrix and solve the matrix system is not compatible with real-time applications. • Angular spectrum projection method [28]: This method's efficiency is quantified as largely faster than the O-mt-TR. At the same time, the obtained target shapes exhibit a very poor resolution and are obtained as a discontinuous collection of discrete points.
Therefore, we conclude that when compared to other available near-field shaping methods, the technique proposed in this work excels in terms of computational efficiency, resolution and suitability. Furthermore, in contrast to most state-of-the-art literature, an accurate quantification of the number of operations, computational time and memory usage were provided.

Conclusions
A novel algorithm was developed to obtain near-field intensity shaping through general antenna array topologies, outperforming traditional techniques. The method leverages the simulated active far-field radiation patterns of the antenna elements. By making use of a spherical Fourier transform and multipole expansion of these radiation patterns, the advocated algorithm performs a far-to-near-field transform to quickly determine the amplitude and phases required to shape the co-polar component of the electromagnetic field produced by general array topologies. Specifically, and as an extension of conventional farfield beam-shaping techniques, the proposed algorithm achieves an accurate shaping over the surface of a sphere with a prescribed radius located in the array's radiative near-field region. The novel method produces sharp patterns, while consuming a very small amount of resources, making it suitable for real-time operation.
As a proof-of-concept, 197-dipoles arrays were exploited to implement near-field focusing for two different intricate shapes. The undesired secondary beams within the target region were significantly suppressed by around 30 dB. Furthermore, intensity-shaping with smaller arrays, comprising 49 dipoles, was also validated using commercial fullwave software. The complex currents required for the shaping and obtained by the novel method were used to excite the full-wave models of the 49 element arrays. A very good agreement between the electric near-fields resulting from both the implemented algorithm and the full-wave calculations was obtained. In an important application scenario, it was shown that the method manages to target individual users who are distributed close to one another in the radiative near-field of the antenna array, as such ensuring the best quality of communication. The method has thus been validated via numerical and fullwave experiments, clearly demonstrating its applicability in (future) B5G and 6G real-time applications and communications scenarios.
Future plans include employing an accurate system impulse response to model applications defined in domains with specific shapes and symmetries. Moreover, an extension to 3D shaping is possible by making use of efficient optimization procedures owing to the algorithm's impressive computation speed.