An Efficient and Accurate Algorithm for Electromagnetic Wave Propagation Modeling Based on Wavelet Transforms

Abstract

In this article, a comprehensive analysis of the approach known as Split-Step Wavelet Parabolic Equation (SSW-PE) in modeling radio wave propagation is presented. The SSW-PE introduces innovations, such as the application of narrow-angle and wideangle approaches, referred to as NAPE and WAPE, respectively. Furthermore, the SSW-PE demonstrates the incorporation of refractivity variations, modifications in terrain modeling for better representation, and considerations of surface boundary conditions. In addition to its innovative aspects, this study aims to provide a complete guide for effectively replicating the algorithm, thereby promoting the advancement of propagation studies using wavelets. The effectiveness and applicability of this approach are validated through comparative studies with well-established solutions, including the Discrete Mixed Fourier Transform (DMFT) version of the Split-Step Parabolic Equation (SSPE) method. Comparisons with measurements from real propagation cases are also conducted. Statistical analysis confirms the innovative potential of the SSW-PE algorithm, which also offers computational efficiency for rapid and consistent simulations. Thus, this article contributes to a comprehensive and innovative analysis, providing tangible resources for the research community interested in expanding this methodology.

Index Terms
Split-Step Wavelet (SSW-PE) Approach; Parabolic equation; Radio Wave Propagation.

I. INTRODUCTION

The primary objective of wireless communication technologies’ evolution is to meet the demand for high-speed transmission and increased communication channel capacity. However, significant challenges are encountered in this evolution due to the use of higher frequencies and broadband communication systems [¹]. When non-homogeneous conditions in environments such as atmospheric variations, irregular terrains, and surface characterization come into play, radio channel characterization becomes even more challenging.

In this context, research in radio propagation models seeks to develop fast, accurate, and refined methods that can capture the complexity of propagation environments [²]-[⁴], including terrain, soil composition, and surface roughness. These models are crucial for predicting signal coverage and characterizing the radio channel in various scenarios, especially in rural and suburban areas where connectivity is essential. Furthermore, they are critical for optimizing antenna positions based on local conditions, predicting the coverage of new communication systems, and assessing the impact of artificial structures such as solar panels or wind turbines on the performance of existing systems.

In the literature, many two-dimensional (2-D) spatial implementations can be found to solve radio propagation problems. Among these, an analytical model is offered as a powerful tool for solving wave propagation problems, with a special focus on truncated Bessel beams in free spaces [⁵]. This approach is not only provided with a deep understanding of the behavior of these beams under various propagation conditions but also presents significant benefits and broad potential applications. By being based on an expansion in Fourier-Bessel series, the model enables precise and efficient theoretical analysis, facilitating the interpretation of experimental results and providing valuable insights across a wide range of fields, including optics and acoustics.

Numerical techniques like the Parabolic Equation (PE) model, Finite Difference (FD), Method of Moments (MoM), Finite Element Method (FEM), and Finite-Difference Time-Domain (FDTD) method have been successfully implemented in both closed and open environments [⁶]. However, accurate modeling of electromagnetic wave propagation in non-homogeneous environments presents significant challenges. Traditional numerical methods [⁷] have limitations related to the required mesh size to ensure accuracy. These methods demand very fine discretizations to achieve accurate results, which can be computationally expensive and impractical for practical scenarios.

The Parabolic Equation Method originated in the 1940s [⁸], initially applied by Leontovich and Fock to address the challenge of vertical stratification in atmospheric propagation. Derivations of the PE from Maxwell’s equations have been extensively discussed in scientific literature [⁷], contributing valuable insights and significant advancements in understanding and modeling various wave propagation problems. However, it is relevant to emphasize that various operator approximations [⁹] can be applied to the Helmholtz wave equation to obtain the PE, and these approximations may result in differences in result accuracy. In this article, we employ a specific formulation [¹⁰] of the standard PE approach, developed to solve two-dimensional problems. This formulation is detailed in (1), where the variables are defined as follows: x represents ground range (longitudinal coordinate), z represents altitude above the Earth’s surface (transverse coordinate), k denotes the free-space wavenumber, n(z) represents the refractive index, and u(x, z) = e^-ikxψ(x, z) signifies the field auxiliary function. The function ψ pertains to the electric or magnetic field in horizontal or vertical polarization, respectively [⁶].

(1) $$\frac{{\mathrm{\partial }}^{2}u(x,z)}{\mathrm{\partial }{z}^2}+2ik\frac{\mathrm{\partial }u(x,z)}{\mathrm{\partial }x}+k^2(n(z{)}^2-1)u(x,z)=0.$$

The Classical Parabolic Equation, as expressed in Equation 1, is transformed into a first-order ordinary differential equation that can be directly solved through a spatial-to-spectral domain Fourier transform. In the context of real-world problems, where the refractive index n can be a function of height and/or range (i.e., n = n(x, z)), this approach is appropriate, as the equation is solved for each small position increment, ∆x, which is chosen sufficiently small to ensure that within any interval ∆x, the refractive index can be considered constant with respect to x. The numerical solution of the NAPE for the SSPE is provided as [⁶], [¹¹], [¹²]:

(2) $$u(x+\mathrm{\Delta }x,z)=\exp \left[\frac{ik\mathrm{\Delta }x}{2}(n^2(x,z)-1)\right]F^{-1}\{\exp \left[\frac{-i{p}^2\mathrm{\Delta }x}{2k}\right]F[u(x,z)]\},$$

where F denotes a Fourier transform and p = k sin θ, where θ is the propagation angle.

Similarly, the numerical solution for the WAPE parabolic equation [¹³ ]-[¹⁵] is defined as follows:

(3) $$\begin{gathered} \begin{array}{}& u(x+\mathrm{\Delta }x,z)=\exp [\frac{ik\mathrm{\Delta }x}{2}(n(x,z)-1)] \\ \times F^{-1}\{\exp \left[\frac{-i{p}^2\mathrm{\Delta }x}{k}{(\sqrt{1-{\left(\frac{p}{k}\right)}^2}+1)}^{-1}\right]F[u(x,z)]\}.\end{array} \end{gathered}$$

In this regard, the Split-Step Wavelet algorithm (SSW) [¹⁶]-[¹⁸] emerges as an innovative alternative in radio channel modeling and characterization. Effectively addressing modeling through the use of a numerical method based on the Parabolic Equations (PE) [⁷], this research aims to make a significant contribution to the characterization of radio channels in realistic scenarios. What makes this research innovative is the incorporation of Wavelet transforms in radio propagation algorithms aimed at analyzing more realistic case studies, in addition to providing channel analysis, with the goal of supporting validation, verification, and confirmation of the technique’s potential.

Wavelets are powerful mathematical tools that enable the decomposition of radio signals into components that can be more efficiently analyzed, contributing to a more accurate channel characterization. This allows for a more efficient signal analysis while retaining relevant information at different scales [¹⁹], [²⁰]. The fast wavelet transform (FWT) is particularly advantageous compared to the fast Fourier transform (FFT) in terms of computational complexity and memory efficiency. This represents a significant advancement in radio channel modeling, enabling a deeper and more effective understanding of propagation phenomena, including their complex interactions with the environment and terrain.

The wavelet-based method also allows for the consideration of multiple factors affecting radio wave propagation, including terrain topography, irregular terrain profiles, and atmospheric refraction. Therefore, the innovation of wavelet-based methods offers a promising approach to modeling radio wave propagation.

This approach ensures fast and efficient simulation times while overcoming the limitations of traditional methods, making it possible to make accurate predictions in both homogeneous and nonhomogeneous scenarios, as well as in realistic environments. Combining the computational efficiency of wavelets with the accuracy of propagation modeling is crucial to meet the demands of modern wireless communication networks.

This study conducts comparative analyses between the proposed algorithm and the Discrete Mixed Fourier Transform (DMFT) version of the Split Step Parabolic Equation (SSPE) method to validate the effectiveness and applicability of the SSW-PE approach. In summary, the use of wavelets in radio channel modeling represents a significant advancement in the topic, enabling the attainment of more accurate and computationally efficient results. With this approach, a valuable solution is provided for the challenges posed by radio wave propagation in more realistic environments, thereby contributing to the development of more effective communication systems and the resolution of complex engineering problems.

A comprehensive description of the SSW-PE algorithm will be provided in Section II, covering the formulation of the parabolic equations technique, the establishment criteria of the model, transmission source modeling, computational domain discretization, environmental constraints, and the crucial pseudocode of the SSW-PE method will be addressed. Subsequently, in Section III, the application of the technique in the engineering analysis of a practical environment with lossy surface parameters and irregular terrain will be demonstrated. The evaluation of the technique’s complexity will be concluded in the same section, followed by final remarks on the research in Section IV.

II. AN IMPROVED SPLIT-STEP WAVELET METHOD

In radio channel modeling and analysis, recent innovation known as the Split-Step Wavelet (SSW) method [²], [¹⁷], [¹⁸] stands out. The SSW-PE approach, adopted in this study, represents a reinterpretation and fresh perspective on the algorithm, characterized by its remarkable capability to handle nonhomogeneous environments. Unlike traditional numerical methods that require very fine discretizations for accuracy, SSW-PE offers an efficient solution that overcomes this limitation. This is achieved through the space-time multiresolution capacity of Wavelets and the adoption of the staircase terrain discretization model. As a result, it becomes possible to consider the various factors outlined in Section I, which impact the modeling of radio wave propagation.

In this section, the aim is to meticulously establish the steps for the development and implementation of the SSW-PE algorithm. Its ability to provide accurate results, even in non-homogeneous scenarios, makes it a valuable tool in developing more effective wireless communication networks and solving complex engineering challenges related to radio propagation. Therefore, it is essential that the discretization process be conducted with great care to avoid potential numerical instabilities. In this context, the steps for implementing the hybrid approach proposed here for the SSW-PE algorithm addressing the advantages of DMFT-SSPE are detailed, as structured in the schematic diagram illustrated in Fig. 1.

Fig. 1
Schematic diagram of the SSW-PE algorithm

For didactic purposes, the elaboration of the sixteen steps of the proposed SSW-PE algorithm is segmented into subsections, namely, sections II-A to II-P, and discussed within each of them.

A. The First Step: Defining and Modeling Source Parameters

In this step, we propose defining a source modeled by the antenna beam pattern, which is characterized by a Gaussian distribution. This modeling approach is widely preferred in most long-range twodimensional radio propagation models. It offers flexibility to adjust beam width and tilt, operating frequency, and polarization, while also providing an approximate representation for parabolic antennas [³], [⁷]. At the beginning of propagation, it is of paramount importance to determine the initial vertical field profile based on the antenna pattern parameters at the initial position (x₀). Factors such as antenna height (x_s) and elevation angle (θ_ELV ) are considered during this source modeling process. In this stage, it is also necessary to specify the range of the observed scenario, encompassing both distance (x_max) and maximum height (z_max) under observation. Providing details related to simulation parameters is crucial, including defining the number of discretizations along the longitudinal and transverse coordinates, represented by N_x and N_z, respectively. Therefore, these values are influenced by the discretization criteria established for the ∆_x and ∆_z steps.

To prevent numerical instabilities, the discretization process must be conducted with caution, taking into account the Nyquist criterion [¹⁰], [²¹]. The numerical stability criterion adopted for the SSW-PE model aims to prevent the amplification of errors and numerical oscillations that could compromise the solution of the PE model. By following this criterion, we ensure stability and obtain computationally reliable and accurate solutions. This is achieved by aligning the discretization step along the transverse coordinates so that the established expression ∆z = π/k sin (θ_max) is satisfied.

Refer to the schematic diagram for details.

B. The Second Step: Modeling the Investigated Terrain

When addressing electromagnetic propagation modeling over irregular terrains using PE, it is possible to simplify the scenario by treating the terrain as a sequence of horizontal steps. This simplification is achieved by discretizing the terrain using the staircase model and accounting for losses. This approach, referred to as the staircase method, is highly efficient and permits the use of larger range steps [¹¹]. The electromagnetic field propagates in segments of constant height, with the suitable conditions applied at the ground surface. When variations in terrain height occur, corner diffraction is omitted, and the field is set to zero on vertical surfaces [⁶]. At this point, we adopt a criterion for relative vertical discretization. In other words, we define N_zⁱ, which relates to the number of discretizations above the ground level in each height segment.

During this stage, it is essential to address the necessary boundary conditions to complete the twodimensional propagation scenario. Thus, we select the appropriate transverse and longitudinal Boundary Conditions (BCs). The region is divided into N_x longitudinal segments and N_z vertical segments. Transverse BCs that need to be satisfied at the Earth’s surface are expressed as [¹⁰]:

(4) $${\alpha }_1\frac{\mathrm{\partial }u(x,z)}{\mathrm{\partial }z}+{\alpha }_{2}u(x,z)=0,$$

where α₁ and α₂ represent constants. For lossy ground, we introduce an impedance BC or the Cauchytype BC, with specific values assigned to α₁ and α₂ depending on the polarization (horizontal or vertical). These values are determined based on the complex relative dielectric constant ϵ, which, in turn, depends on the relative permittivity (ϵ_r) and the conductivity (σ_g) of the lossy ground surface at each range step [⁶].

Refer to the schematic diagram for details.

C. The Third Step: Defining the Discretized Mesh of the Environment

The domain is discretized considering the parameters N_x and N_z. Consequently, at the initial propagation position x₀, the field to be propagated is known and must be determined through analytical modeling. Subsequently, this field is propagated along the positions of interest p_x through an iterative process until reaching the maximum distance x_max and the maximum height under observation z_max.

To achieve this goal, the domain is discretized along the x and z axes, as defined below:

(5) $$\begin{array}{rl}{x}_{p_x}& =x_0+p_x\mathrm{\Delta }x,for{p}_x\in \{0,\dots ,N_x-1\},\\ z_{p_z}& =p_z\mathrm{\Delta }z,for{p}_z\in \{0,\dots ,N_z-1\},\end{array}$$

where N_x = (x_max - x_s)/∆x and N_z = z_max/∆z represent the numbers of sampling points along the x and z axes, respectively. Here, ∆x and ∆z represent the step sizes along the x and z axes, as mentioned in the second stage.

Once the domain discretization is established, the discretized version of the field profile u(x, z) along z is stored and denoted as u_x[p_z]. In order to ensure the technique’s accuracy and avoid numerical instabilities, the discretization process must be conducted carefully, taking into account the Nyquist criterion, which was revisited in the second stage. The adopted numerical stability criterion for the SSW-PE model aims to prevent the amplification of errors and numerical oscillations that could compromise the solution of the PE model. By adhering to this criterion, we ensure stability and obtain computationally reliable and accurate solutions.

Refer to the schematic diagram for details.

D. The Fourth Step: Defining Wavelet Parameters and Configuring the Scattering Matrix

In this stage, the necessary information for performing the wavelet transformations is introduced, which are selected based on criteria such as compression, regularity, symmetry, and family choice [²²].

There are two possible approaches to define SM as detailed in [¹⁷] and [¹⁸]. In [¹⁷], the matrix is defined with dimensions of (N_z, N_z), referred to as mSSW, containing all propagation elements. Meanwhile, in [¹⁸], the dimensions are (2^L, N_z), referred to as lSSW, and include only non-redundant elements. This latter approach utilizes Wavelet properties to add versatility to the method and reduce storage requirements. In practice, memory usage for storing SM is significantly reduced [¹⁸] due to the small size of the Wavelet support (N_z^s) for the observed domain (N_z^s ≪ N_z). In this study, we adopt the lSSW approach in conjunction with computational implementation details used by the DMFT-SSPE numerical algorithm, aiming to mitigate numerical instability issues, minimize memory usage, and improve computational efficiency.

Algorithm I: Computing the scattering matrix SM.

Input:L - Maximum level of decomposition Wavelet

Output:SM - scattering matrix

1 forl^′ ∈ [1, L] do

2$A_{(l^{\mathrm{\prime }},p_z)}^U⟵$ Applying FWT to a null wavelet vector sized as compact support;

3$A_{(l^{\mathrm{\prime }},p_z)}^U⟵$ set to unit amplitude at central position;

4${\chi }_{(l,0)}[p_z]⟵IFWT{A}_{(l^{\mathrm{\prime }},p_z)}^U$;

5${\chi }_{(l^{\mathrm{\prime }},\mathrm{\Delta }x)}[p_z]⟵$ receives the DMFT-SSPE propagation of ${\chi }_{(l,0)}[p_z]$ by Δx;

6$q_{max}⟵2^{L-l^{\mathrm{\prime }}}$ Number of operators required for level;

7 forp_t ∈ [1, q_max] do

8χ_(l′,∆x)[p_z + p_t] ←- by translating χ_(l′,∆x)[p_z] over p_t points;

9SM_l^′,pt ←- application of FWT and compression V_p on χ_(l′,∆x)[p_z + p_t];

The calculation of SM should be carried out following the proposed systematic procedure, which consists of the following steps, as outlined in the algorithm presented in Algorithm I:

Firstly, The library $A_{(l^{\mathrm{\prime }},p_z)}^U$ of Wavelet coefficients is obtained. For each resolution level, i.e., for all values of l^′ ∈ [1, L], the Wavelet coefficients are determined by defining a null auxiliary vector with a size equivalent to the size of the compact support of the predefined Wavelet family and level. This vector is then subjected to Wavelet transformation, resulting in a sparse vector. Subsequently, the central position of this vector at each Wavelet level is set to unity. After this step, the signal is restored through inverse transformation, meaning that the Wavelet is determined at position χ_(l′,0) along p_z ∈ [0, N_z - 1] by inversely transforming a set of null Wavelet coefficients at all positions of $A_{(l^{\mathrm{\prime }},p_z)}^U$, except at the central position ($A_{(l^{\mathrm{\prime }},0)}^U$), where it is defined with unit amplitude before the application of the Inverse FWT (IFWT).

Next, the library of Wavelet coefficients χ_(l′,0)[p_z] obtained should be propagated over ∆x for each Wavelet level. The maximum level of decomposition Wavelet q_max is responsible for defining the number of signal shifts at each level using the DMFT-SSPE method.

Subsequently, the Wavelets are shifted as a combination of matrices restricted by the scattering matrix SM in order to reduce the number of required Fourier transforms. Thus, each Wavelet propagated at level l^′ is translated from p_t ∈ [0, 2^L-l′ ] along p_z to obtain χ_(l′,∆x)[p_z + p_t] for all values of p_z + p_t ∈ [0, N_z - 1].

Finally, the Wavelets χ_(l′,∆x)[p_z + p_t] are decomposed and compressed with the threshold V_p to obtain SM_l′,pt.

To achieve this objective, the advantages of Wavelet properties [¹⁹], [²²] are exploited, which include the use of temporal and scale domain shifts, scale and translation invariance, energy preservation, and multiresolution decomposition. In obtaining $A_{(l^{\mathrm{\prime }},p_z)}^U$, the compact support of Wavelets is utilized, resulting in a vector dimension significantly smaller than the discretization in p_z. Furthermore, the property leading to a highly sparse matrix containing only significant coefficients is leveraged, thus reducing information storage

Refer to the schematic diagram for details.

E. The Fifth Step: Obtaining the Refractivity Profile

The refractivity profile is obtained by considering the atmospheric conditions of the surface duct and can be modeled throughout the domain in N_z. The refractive index n is related to the dielectric constant ϵ_r of the medium, such that $n=\sqrt{{\mu }_r{ϵ}_r}$, where the (relative) permeability µ_r is defined as one for air. The refractive index of air at the Earth’s surface is approximately n ≈ 1.0003. For convenience, the refractivity profile N , defined as N = (n - 1) × 10⁶ in N-units (dimensionless). The refractive index of air decreases with altitude as the air becomes less dense, exhibiting an approximately linear dependence in the lower part of the atmosphere.

In a standard atmosphere [²³], the refractivity gradient is dn/dh = -40N/km for long-distance observations. The standard atmosphere is characterized by a decrease in temperature, pressure, and humidity with altitude. Although atmospheric pressure generally decreases with altitude, anomalies in the temperature profile can occur, including rapid decreases and increases at specific altitude ranges. As a result, layers with extremely negative refractivity gradients can occur. Another effect leading to extremely high refractivity gradients is caused by the rapid decrease in humidity above large bodies of water.

Due to these phenomena, researchers have conducted numerous studies in search of improved modeling methods. We encourage readers to explore alternative research approaches within the scope of the proposed methodology.

Refer to the schematic diagram for details.

F. The Sixth Step: Obtaining the Attenuation Filter

Due to the nature of the propagation problem, which encompasses a vertically open region, it is necessary to impose boundary conditions at the upper termination of the computational numerical domain. The Sommerfeld radiation condition, where (u_x[p_z])_pz→∞ = 0, must be satisfied, implying that radio waves are completely absorbed upon reaching the upper boundary, and non-physical reflections from this boundary should not interfere with the computational domain [²⁴].

In order to avoid undesirable reflections that may impact numerical results, it is necessary to gradually reduce the field at the upper boundary, avoiding abrupt truncation of the field along the vertical coordinate z. This gradual reduction can be achieved through various approaches, such as the use of absorption regions, the application of windowing functions (such as the Tukey window or the Hanning window), or the introduction of an imaginary part in the refractive index to cause attenuation of radio waves. In this work, we implement and calculate a Hanning window according to the method proposed by [²⁴]:

(6) $$\begin{array}{c}J[p_z]=\{\begin{array}{ll}\frac{1+cos\left(\pi \left[\frac{p_z}{N_z^H}\right]\right)}{2},& p_z\in [N_z-N_z^H,N_z]\\ 1,& p_z\in [0,N_z-N_z^H],\end{array}\end{array}$$

where N_z^H represents the total width of the absorption layer. Mathematically, this windowing function ensures that J[p_z ≤ (N_z - N_z^H)] = 1 and J(N_z) = 0. An example is illustrated in Fig. 2.

Fig. 2
Hanning Windows.

Refer to the schematic diagram for details.

G. The Seventh Step: Initiating the Repetition Cycle for the Auxiliary Field Function

At this stage, the repetition loop is initiated, ranging from p_x = x₀ to N_x. Here, the auxiliary field function u_x[p_zⁱ ] is determined based on the auxiliary field function, with u₀ assigned when p_x = 0, or u_x for other values of p_x. This allocation takes into account the conditions related to field displacement to the modeled position, applying the Staircase approach to represent the terrain at p_x. Simultaneously, the local refractive index n_px = n[N_z - N_zⁱ: N_z] is pre-allocated to match the corresponding position in the terrain profile.

Refer to the schematic diagram for details.

H. The Eighth Step: Applying Refraction Parameters

(7) $$\begin{array}{c}v[p_z^i]=\{\begin{array}{ll}u[p_z^i]\exp (ik((n_{p_x}{)}^2-1)\frac{\mathrm{\Delta }x}{2}),& \text{for NAPE,}\\ u[p_z^i]\exp (ik(n_{p_x}-1){\mathrm{\Delta }}_x),& \text{for WAPE.}\end{array}\end{array}$$

Refer to the schematic diagram for details.

I. The Ninth Step: Applying the Discrete Approach for the coefficients

The discretization based on [²⁵], with a weighted factor of 0.5 used for the first and last terms (i.e., p_zⁱ = 0, p_zⁱ = N_zⁱ), to obtain coefficients C₁, C₂, r, and D.

The following quantities are required to complete the definition of the transform:

(8) $$C_1[p_z^i]=D\sum _{p_z^i=0}^{N_z^i}v[p_z^i]r^{p_z^i},$$

(9) $$C_2[p_z^i]=D\sum _{p_z^i=0}^{N_z^i}v[N_z^i-p_z^i](-r{)}^{N_z^i},$$

where,

(10) $$D=\frac{2(1-r^2)}{(1+r^2)(1-r^{2(N_z^i)})},$$

and r, -r^-1 are the roots of the quadratic equation:

(11) $$r^2+2r\alpha \mathrm{\Delta }z-1=0,$$

where r depends on the polarization of the transmitter:

(12) $$\begin{array}{c}r=\{\begin{array}{ll}\sqrt{1+(\alpha \mathrm{\Delta }z{)}^2}-\alpha \mathrm{\Delta }z,& \text{for vertical,}\\ -\sqrt{1+(\alpha \mathrm{\Delta }z{)}^2}-\alpha \mathrm{\Delta }z,& \text{for horizontal.}\end{array}\end{array}$$

Refer to the schematic diagram for details.

J. The Tenth Step: Constructing the auxiliary field

The adopted method is based on the implementation proposal of DMFT as presented in [²⁶], where the vector of the function w[p_zⁱ ] is defined as follows:

(13) $$w[p_z^i]\equiv \frac{v[p_z^i+1]-v[p_z^i-1]}{2\mathrm{\Delta }z}+\alpha v[p_z^i],$$

where p_zⁱ belongs to the interval [1, N_zⁱ - 1].

Refer to the schematic diagram for details.

K. The Eleventh Step: Wavelet Propagation

In this context, the modified function w[p_zⁱ ] undergoes propagation using the scattering matrix. The calculation of w[p_zⁱ ] follows a careful procedure as indicated in the algorithm in Fig. II.

First, we apply the Fast Wavelet Transform and compress the data V_s into the auxiliary function w[p_zⁱ ]. Next, we proceed to dilate w[p_zⁱ]FWT for each decomposition level l^′ as required for the respective level. Propagation in ∆x is performed through the convolution of the relevant elements in SM_l,qmax[l′]. Finally, we apply the Inverse Fast Wavelet Transform to the convolution result, allowing us to establish w[p_zⁱ , ∆x].

Refer to the schematic diagram for details.

Algorithm II: propagated field without effects of the medium.

input :L - Maximum level of decomposition Wavelet,

SM - Scattering Matrix,

w[p_zⁱ ] - field function,

q_max - number of required operators by level

output:w[p_zⁱ , ∆x] - propagated field without effects of the medium

1W [p_zⁱ ] ← CWw[p_zⁱ ] wavelet decomposition of the field and compression (V_s);

2 foreachl^′ ∈ [1, L] do

3nel ← q_max[l^′]× number of non-zero coefficient of W [p_zⁱ ];

4U_dilated ← array of zeros of size nel;

5U_dilated[(1 : q_max[l^′] : nel)] ← W [p_zⁱ ] coefficients Wavelets dilated in q positions;

6W [p_zⁱ , ∆x][l^′] ← W [p_zⁱ , ∆x][l^′] + (SM_{l,qmax [l^′]} ∗ U_dilated) by convolution ;

7w[p_zⁱ , ∆x] ← inverse transformation Wavelet of W [p_zⁱ , ∆x];

L. The Twelfth Step: Multiplying by the Spatial-p Operator

The multiplicative process by the associated spatial-p operator is applied based on the chosen propagation solution. It is performed as follows:

(14) $$\begin{array}{c}wp[p_z^i,\mathrm{\Delta }x]=\{\begin{array}{ll}w[p_z^i,\mathrm{\Delta }x]\exp \left(i{p}^2\frac{\mathrm{\Delta }x}{2k}\right),& \text{for NAPE,}\\ w[p_z^i,\mathrm{\Delta }x]\exp (-\frac{i{p}^2\mathrm{\Delta }x}{k}{(\sqrt{1-\frac{p^2}{k^2}}+1)}^{-1}),& \text{for WAPE,}\end{array}\end{array}$$

where p is a vector of size N_zⁱ, and the vector elements are generated with equal increments, starting from p_max/N_zⁱ and ending at p_max, with p_max = k sin θ_max.

Refer to the schematic diagram for details.

M. The Thirteenth Step: Propagating Coefficients

The process consists of performing the following procedure:

(15) $$C_1[p_z^i,\mathrm{\Delta }x]=C_1[p_z^i]\exp \left(\frac{i\mathrm{\Delta }x}{2k}{\left(\frac{\log r}{\mathrm{\Delta }z}\right)}^2\right),$$

(16) $$C_2[p_z^i,\mathrm{\Delta }x]=C_2[p_z^i]\exp \left(\frac{i\mathrm{\Delta }x}{2k}{\left(\frac{\log (-r)}{\mathrm{\Delta }z}\right)}^2\right).$$

Refer to the schematic diagram for details.

N. The Fourteenth Step: Obtaining the Particular Solution

The solution is a combination of the particular solution of the inhomogeneous equation and an appropriate linear combination of the general solutions of the homogeneous equation:

(17) $$v[p_z^i,\mathrm{\Delta }x]=v_p[p_z^i]+A{r}^{p_z^i}+B(-r{)}^{(N_z^i-p_z^i)}.$$

The particular solution v_p[p_zⁱ ] is found by factoring the equation of differences (13) into upper and lower triangular terms and using the recursions:

(18) $$y[p_z^i]-ry[p_z^i-1]=2\mathrm{\Delta }zwp[p_z^i,\mathrm{\Delta }x],$$

(19) $$v_p[p_z^i+1]+\frac{1}{r}{v}_p[p_z^i]=y[p_z^i],$$

for p_zⁱ = 1, ..., N_zⁱ - 1, with y[0] and v_p[N_zⁱ] being arbitrary. For Eq. 18, the solution is retro-resolved. Direct substitution verifies that these equations are equivalent to Eq. 13. The coefficients A and B are defined as:

(20) $$A=C_1[p_z^i,\mathrm{\Delta }x]D\sum _{p_z^i=0}^{N_z^i}{v}_p[p_z^i]r^{p_z^i},$$

(21) $$B=C_2[p_z^i,\mathrm{\Delta }x]D\sum _{p_z^i=0}^{N_z^i}{v}_p[N_z^i-p_z^i](-r{)}^{N_z^i}.$$

Refer to the schematic diagram for details.

O. The Fifteenth Step: The Final Calculation Phase in the Repetition Cycle

In this phase, the refractivity parameters and the Hanning Window are applied by the spatial operator associated with the refractive index, according to the adopted solution for propagation, as indicated below:

(22) $$\begin{array}{c}u[p_z^i,\mathrm{\Delta }x]=\{\begin{array}{ll}v[p_z^i,\mathrm{\Delta }x]\exp (ik(n^2-1)\frac{\mathrm{\Delta }x}{2})J[p_z^i],& \text{for NAPE,}\\ v[p_z^i,\mathrm{\Delta }x]\exp (ik(n-1){\mathrm{\Delta }}_x)J[p_z^i],& \text{for WAPE.}\end{array}\end{array}$$

Next, the value of u[p_zⁱ , ∆x] is assigned to both u_x and the auxiliary vector that accumulates the propagation positions, u[p_x, p_z].

Once this process is completed, it is considered that the information related to the propagation in the p_x step is finalized, and the index p_x is incremented by one position, that is, p_x + 1 = ∆x + ∆x. Refer to the schematic diagram for details.

P. The Sixteenth Step: The Decision-Making Process Related to Propagation

At this point, it is evaluated whether p_x has reached the final position x_max of the environment under analysis. If not, the repetition cycle is automatically restarted, returning to the seventh stage. If p_x has reached the final position, the auxiliary field vector u[p_x, p_z] is concluded and forwarded. This vector can then be processed for the purpose of generating reports, including creating graphs, analyzing path losses, and calculating power.

Refer to the schematic diagram for details.

III. CASE STUDIES

Case studies are conducted to address the formulation and the presented SSW-PE algorithm in Section II.

A. Practical Case: Jerslev Measurement Campaign

In this practical case, the measurement campaign [²⁷], [²⁸] took place in Northem Jutland, near Aalborg, Denmark. The environment, named Jerslev, under analysis, exhibits an irregular profile with rural areas, trees, and some buildings, spanning approximately 5.5 km. Terrain data were acquired from digital maps with a 50 m resolution, as illustrated in Fig. 3. During the measurement campaign, a half-wave dipole with a height of 10.4 m and vertical polarization was employed as the transmitter, radiating 10 W at a signal frequency of 144 MHz. A receiving station, positioned along the terrain at a height of 2.4 m, recorded the electric field intensity in the sample space.

Fig. 3
Comparison of Measured and Estimated Path Loss in Different Techniques at Jerslev, Northen Jutland, Denmark.

Standard electrical properties for the soil characterization were used to obtain numerical method results, with ϵ_r = 15 and σ = 0.012 S/m. Following discretization criteria, ∆x = 50 m and ∆z = 1 m were employed. A vertical domain limit was set at z_max = 200 m. The Fig. 3 illustrates the measured path loss, compared among the MoM, FDTD, SSPE-NAPE, SSPE-WAPE, SSW-PE-NAPE, and SSPW-PE-WAPE techniques along the terrain profile.

Table I presents error statistics for the case, with measurements serving as the reference for calculating mean absolute error and standard deviation. Furthermore, the dispersion between the numerical techniques and the measurements is illustrated in Fig. 4. Data analysis confirms that the numerical techniques yield similar results, and the SSPW-PE algorithm provides a similar solution in both scenarios, whether using the NAPE or WAPE model. In addition to the statistical comparison, a computational efficiency assessment is conducted, revealing the processing time required to obtain results with each technique. This strengthens and validates the effectiveness of the SSW-PE technique, demonstrating its ability to achieve good convergence of results and remarkable computational efficiency compared to other numerical approaches.

TABLE I
Statistical Parameters for the Jerslev Case Study

Method	Mean value error[dB]	Standard deviation [dB]	Simulation time
FDTD	4.1554	5.0334	50.3 hours
MoM	3.0226	3.7738	6.9 hours
NAPE	4.0079	4.5595	31.72 seconds
SSW-PE-NAPE	8.5287	10.1518	0.57 seconds
SSW-PE-WAPE	7.3020	8.8188	0.89 seconds
WAPE	4.0785	4.6369	31.72 seconds

Fig. 4
Scatter plot comparing Path Loss Estimates between Different Techniques and Measurements at Jerslev, Northern Jutland, Denmark.

B. Practical Case: Campaign

In this specific scenario, the measurement campaign [²⁹] was conducted within the metro-rail transportation system in Belo Horizonte, Minas Gerais, Brazil (BH/MG/BRA). This system extends for approximately 28 km and encompasses 19 passenger boarding and disembarking stations. In 2019, the Brazilian Urban Trains Company (CBTU) conducted a radio propagation study to assess the coverage of this new system along the railway.

During the field survey phase, eight potential locations for the installation of Base Radio Stations (BRS) were identified. Subsequently, a drive test was executed in which a transmitter was placed at the predetermined BRS locations, while a reception system was installed on a passenger train. This system recorded georeferenced data regarding received signal power as the train traversed the railway. The transmission system utilized a vertical monopole antenna with a gain of 8.15 dBi, connected to a bidirectional repeater for digital radio with a maximum transmission power capacity of 50 W.

The measurement campaign methodology was designed to operate the transmission system at a frequency of 163.94625 MHz with a power output of 20 W. The reception system recorded the received signal level every 10 seconds at a height of 2 meters as the train traveled along the railway. The train made stops at each station along the route to facilitate passenger boarding and disembarking. The collected data is available in [²⁹].

Elevation data for the terrain relevant to the analyzed links was obtained from data provided by Embrapa Monitoring by Satellite [³⁰], through the Brazil in Relief project. The digital elevation models have a spatial resolution of 90 meters and use the WGS-84 geographic coordinate system. This data was processed using geoprocessing software, and the terrain was discretized with a horizontal spacing of 12 meters.

For the purpose of analysis in this article, we present the link between the Radio Base Station (ERB), located at the transmission station called Minas Shopping, and the distribution of the electric field, as well as the attenuation along the path directed to São Gabriel Station, as illustrated in Fig. 5 and Fig.6, respectively. In the case of losses along the path, a comparison between the DMFT-SSPE-NAPE, DMFT-SSPE-WAPE, SSW-PE-NAPE, and SSPW-PE-WAPE techniques is provided. To obtain results through numerical methods, standard soil electrical properties with a relative permittivity (ɛ_r) value of 15 and conductivity (σ) of 0.012 S/m were used for characterization. According to discretization criteria, the values of ∆z were aligned with the Nyquist criterion, as described in Section II-A. A vertical domain limit (z_max) was defined at 2000 meters.

Fig. 5
Normalized electric field distribution along the path initiated at Minas Shopping Metro Station towards São Gabriel Station, BH/MG/BRA.

Upon data analysis, it was observed that the results obtained by the numerical techniques DMFT- SSPE and SSW-PE are similar. As expected, both the DMFT-SSPE and SSW-PE models produce nearly identical results when considering both NAPE and WAPE solutions. The overlap of data obtained through narrow-angle and wide-angle approaches for both numerical solution techniques mentioned earlier is illustrated in Fig. 6. Subsequently, information regarding the processing time of the techniques is presented in Table II, highlighting the computational efficiency of the SSW-PE algorithm. The investigation is concluded with the analysis of the vertical propagation profile at the receiver’s location, namely São Gabriel Station, at a distance of 1015 meters from the transmitter, where a received power of -76 dBm was measured during the measurement campaign [²⁹]. The received power for this vertical profile at the mentioned distance is illustrated in Fig. 7, where an initial oscillation of approximately 8 dB is observed in the numerical techniques.

TABLE II
Comparison of processing time among different numerical techniques for the directed path from Minas Shopping Metro Station towards São Gabriel, BH/MG/BRA.

Method	Simulation time
DMFT-SSPE-NAPE	10.56 seconds
DMFT-SSPE-WAPE	24.77 seconds
SSW-PE-NAPE	1.71 seconds
SSW-PE-WAPE	2.10 seconds

Fig. 6
Comparison of path loss along the route among different numerical techniques for a horizontal profile starting at Minas Shopping Metro Station towards São Gabriel Station, BH/MG/BRA.

Fig. 7
Comparison of received power among different numerical techniques for a vertical profile located at São Gabriel Metro Station, BH/MG/BRA.

The specific position chosen for comparison with the measurements was determined based on the nature of the reference work [²⁹], where the measurement campaign was conducted at specific transmission points determined by the transmitting stations of the urban train service in the city of Belo Horizonte, MG, Brazil. Given the transmission profile established between the investigated locations, it was possible to compare with the point arranged in the analysis of the vertical transmission profile.

IV. CONCLUSION

In this article, the effectiveness of the SSW-PE algorithm is highlighted as an efficient and accurate alternative for electromagnetic wave propagation predictions. A comprehensive guide to the SSW-PE has been presented, comprising a series of detailed steps that enable and encourage the scientific community to engage in investigating the advantages of using Wavelet transforms compared to the established Fourier-based models. Realistic studies conducted here allowed the observation of results obtained by the SSW-PE in comparison to other numerical techniques and measurement campaigns. In case studies conducted in Jerslev, Denmark, and the metro-rail transport system located in Belo Horizonte, Minas Gerais, Brazil, it was confirmed that SSW-PE successfully reproduced path losses in lossy terrains, irrespective of the propagation model used (NAPE or WAPE). Furthermore, it was demonstrated that the algorithm is characterized by its computational efficiency, requiring significantly less processing time compared to other numerical approaches.

In summary, the utility of the SSW-PE as a valuable tool for wave propagation modeling is emphasized, with the potential for application in various real-world scenarios. The SSW-PE offers a combination of accuracy and efficiency, making it suitable for a wide range of applications in wireless communications, radar systems, and other related fields.

ACKNOWLEDGMENTS

This work was carried out with the support of CNPq (under grant number 311521/2020-7) and FAPEMIG (under grant APQ-02286-23).

Our sincere gratitude is expressed to Rodrigo Oliveira, a Master of Electrical Engineering from the Graduate Program in Electrical Engineering at the Federal University of Minas Gerais, for the generous contribution of data from the metro-rail transportation system located in Belo Horizonte, Minas Gerais, Brazil, which played a crucial role in the completion of this study.

REFERENCES

• ^[1]
  T. S. Rappaport, Y. Xing, G. R. MacCartney, A. F. Molisch, E. Mellios, and J. Zhang, “Overview of millimeter wave communications for fifth-generation (5g) wireless networks-with a focus on propagation models,” IEEE Transactions on Antennas and Propagation, vol. 65, no. 12, pp. 6213-6230, 2017.
• ^[2]
  H. Zhou, A. Chabory, and R. Douvenot, “A Fast Wavelet-to-Wavelet Propagation Method for the Simulation of Long- Range Propagation in Low Troposphere,” IEEE Transactions on Antennas and Propagation, vol. 70, no. 3, pp. 2137-2148, 2022.
• ^[3]
  D. Parada, C. G. Rego, D. Guevara, A. Navarro, G. L. Ramos, and R. Oliveira, “A modified radiopropagation multipath model for constant refractivity gradient profiles,” Journal of Microwaves, Optoelectronics and Electromagnetic Applications, vol. 22, no. 2, p. 298-312, 2023.
• ^[4]
  G. Dockery, “Modeling electromagnetic wave propagation in the troposphere using the parabolic equation,” IEEE Transactions on Antennas and Propagation, vol. 36, no. 10, pp. 1464-1470, 1988.
• ^[5]
  C. Dartora and K. Nobrega, “Study of gaussian and bessel beam propagation using a new analytic approach,” Optics Communications, vol. 285, no. 5, pp. 510-516, 2012.
• ^[6]
  G. Apaydin and L. Sevgi, Radio wave propagation and parabolic equation modeling John Wiley & Sons, Ltd, 2017.
• ^[7]
  P. Zhang, L. Bai, Z. Wu, and L. Guo, “Applying the parabolic equation to tropospheric groundwave propagation: A review of recent achievements and significant milestones.” IEEE Antennas and Propagation Magazine, vol. 58, no. 3, pp. 31-44, 2016.
• ^[8]
  M. Leontovich and V. Fzck", “Solution of propagation of electromagnetic waves along the earth’s surface by the method of parabolic equations,” Journal of Physics USSR vol. 10 pp. 13-23, 1946.
• ^[9]
  D. J. Thomson and N. R. Chapman, “A wide-angle split-step algorithm for the parabolic equation,” The Journal of the Acoustical Society of America, vol. 74, no. 6, pp. 1848-1854, 1983.
• ^[10]
  O. Ozgun, G. Apaydin, M. Kuzuoglu, and L. Sevgi, “Petool: Matlab-based one-way and two-way split-step parabolic equation tool for radiowave propagation over variable terrain,” Computer Physics Communications, vol. 182, no. 12, pp. 2638-2654, 2011.
• ^[11]
  M. Levy, Parabolic Equation Methods for Electromagnetic Wave Propagation IET, 2000, vol. 45.
• ^[12]
  J. R. Kuttler and G. D. Dockery, “Theoretical description of the parabolic approximation/fourier split-step method of representing electromagnetic propagation in the troposphere,” Radio Science, vol. 26, no. 02, pp. 381-393, 1991.
• ^[13]
  J. F. Claerbout, “Fundamentals of geophysical data processing,” International Series in the Earth and Planetary Sciences, 1985.
• ^[14]
  M. D. Feit and J. A. Fleck, “Light propagation in graded-index optical fibers,” Applied Optics, vol. 17, no. 24, pp. 3990-3998, 1978.
• ^[15]
  M. D. Collins and R. B. Evans, “A two-way parabolic equation for acoustic backscattering in the ocean,” The Journal of the Acoustical Society of America, vol. 91, no. 3, pp. 1357-1368, 1992.
• ^[16]
  A. Iqbal and V. Jeoti, “A split step wavelet method for radiowave propagation modelling in tropospheric ducts,” in 2011 IEEE International RF & Microwave Conference, 2011.
• ^[17]
  H. Zhou, R. Douvenot, and A. Chabory, “Modeling the long-range wave propagation by a split-step wavelet method,” Journal of Computational Physics, vol. 402, p. 109042, 2020.
• ^[18]
  T. Bonnafont, R. Douvenot, and A. Chabory, “A local split-step wavelet method for the long range propagation simulation in 2d,” Radio Science, vol. 56, no. 2, pp. 1-11, 2021.
• ^[19]
  S. Mallat, “A theory for multiresolution signal decomposition: the wavelet representation,” IEEE Transactions on Pattern Analysis and Machine Intelligence, vol. 11, no. 7, pp. 674-693, 1989.
• ^[20]
  I. Daubechies, “Orthonormal bases of compactly supported wavelets,” Communications on Pure and Applied Mathematics, vol. 41, no. 7, pp. 909-996, 1988.
• ^[21]
  O. Ozgun, M. Kuzuoglu, G. Apaydin, and L. Sevgi, “Two-way split-step parabolic equation algorithm for tropospheric propagation: Tests and comparisons,” in 2010 10th Mediterranean Microwave Symposium, pp. 14-17, 2010.
• ^[22]
  S. Mallat, A Wavelet Tour of Signal Processing, Third Edition: The Sparse Way, 3rd ed. USA: Academic Press, Inc., 2008.
• ^[23]
  I. C. Union, ITU-R P.453-14 The radio refractive index: its formula and refractivity data, ITU Recommendation”, 2019.
• ^[24]
  P. Zhang, L. Bai, Z. Wu, and F. Li, “Effect of window function on absorbing layers top boundary in parabolic equation,” in Asia-Pacific Conference on Antennas and Propagation, pp. 849-852, 2014.
• ^[25]
  D. Dockery and J. Kuttler, “An improved impedance-boundary algorithm for fourier split-step solutions of the parabolic wave equation,” IEEE Transactions on Antennas and Propagation, vol. 44, no. 12, pp. 1592-1599, 1996.
• ^[26]
  J. R. Kuttler and R. Janaswamy, “Improved fourier transform methods for solving the parabolic wave equation,” Radio Science, vol. 37, no. 2, pp. 1-11, 2002.
• ^[27]
  J. Hviid, J. Andersen, J. Toftgard, and J. Bojer, “Terrain-based propagation model for rural area-an integral equation approach,” IEEE Transactions on Antennas and Propagation, vol. 43, no. 1, pp. 41-46, 1995.
• ^[28]
  C. Garcia Batista and C. Gonçalves do Rego, “A high-order unconditionally stable fdtd-based propagation method,” IEEE Antennas and Wireless Propagation Letters, vol. 12, pp. 809-812, 2013.
• ^[29]
  R. Oliveira, G. L. Ramos, and C. G. do Rego, “Predição de cobertura radioelétrica usando métodos de equação parabólica em ambiente metroviário urbano na faixa de vhf,” in Programa de Pós-Graduação em Engenharia Elétrica - PPGEE, 2021.
• ^[30]
  E. Territorial, Satélites de Monitoramento, 2023.

Publication Dates

• Publication in this collection
  15 Apr 2024
• Date of issue
  2024

History

• Received
  19 Oct 2023
• Reviewed
  23 Nov 2023
• Accepted
  21 Feb 2024