1 Introduction

An important part of geophysical studies is to conclude the inner earth from data collected at or near the surface of the earth. The measured data is somehow connected to one or more interesting properties of the earth. These physical properties can be adequately estimated by solving an inverse problem. After modelling the subsurface volume under investigation, an inverse problem can be resolved to find the probable values of the model parameters. In geophysical inversion, the objective is to find an earth model defined by a set of physical values that matches observational data. Data acquisition techniques are undergoing permanent improvements over the years, hence, more cutting-edge data processing methods are required. The objective of this study was to develop new data processing algorithms by applying inverse theory to Fourier transformation for geophysical applications. In signal processing, a change of data from the time domain to frequency domain is a common practice. Commonly, the Discrete Fourier Transformation (DFT) algorithm is used to establish the discrete frequency components of a distinctly sampled time-domain dataset. However, the noise sensitivity of the processing method must be taken into account as the noise recorded in the time domain is adequately transformed into the frequency domain. In order to provide significant reduction of the impact of data noise and outliers in data processing, a more robust method known as the Steiner Iteratively Reweighted Least Square Fourier Transformation (S-IRLS-FT) was developed by Dobróka et al. (2015). The algorithm was evaluated on several levels including an application on magnetic data set (Dobróka et al. 2017) to prove its noise reduction capability. By combining the IRLS algorithm with Cauchy-Steiner weights the Fourier transform was handled as a robust inverse problem and the Fourier spectrum discretized utilizing series expansion. These series expansion based inversion methods have been successfully applied in processing gravity data (Dobróka and Völgyesi 2010), borehole geophysical data (Szabó 2004; Dobróka et al. 2016) as well as Induced polarization data (Turai et al. 2010).

2 Theoretical basis of the IRLS‑FT method

The developed algorithm uses series expansion based discretization of the Fourier spectrum with Legendre polynomials as basis functions of discretization, and the solution of an inverse problem provides the estimated values of expansion coefficients. The robustified process uses a variant of the IRLS method, which is relied on the Steiner weight instead of the Cauchy one to facilitate the internal iterative recalculation of weights.

For a quicker and easier solution of an inverse problem, we normally apply the simple least-squares method, which is very effective when data noise follows Gaussian distribution. In situations where the data errors are not consistent and sparsely distributed such as outliers, the calculated model value may be highly questionable, which creates a limiting factor to the application of the least-squares method since geophysical datasets often contain outliers (Nuamah and Dobroka 2019). Several theoretically accepted methods have been formulated over the years to adequately process data outliers in other to achieve statistical robustness. The Least Absolute Deviation (LAD) as a robust method minimizes the L1 norm of the deviation vector characterizing the misfit function between the observed and predicted data. Continual practical affirm investigations that inversion with minimization of the L1 norm provides good estimates when fewer large errors contaminate the data (Dobroka et al. 2017). An easy approach is to use the Cauchy criterion which adopts a random-distributed data noise. To properly account for the variations in data noise during inversion, each data must contribute to the solution based on its error margins and this is achieved by weighting the data. By applying Cauchy weights, Cauchy inversion can be used as a robust optimization method (Amundsen 1991). Although a very efficient procedure is obtained by integrating the IRLS algorithm with Cauchy weights, it is a little challenging since the scale parameter of the weights has to be a-priori known. This challenge was resolved by the Most Frequent Value (MFV) method developed by Steiner (1988,1997) which derives the scale parameters from the real statistics of the data set. Dobróka et al. 1991 first emphasized the possibility of inserting calculated MFV-weights based on Steiner’s method into an IRLS algorithm to produce a very efficient and robust method. A successful application of the MFV method in processing borehole geophysical data was reported by Szűcs et al. (2006). The so-called Cauchy-Steiner weights were also useful in robust tomography algorithms by Szegedi and Dobróka (2014).

2.1 Legendre polynomials as basis functions

Legendre polynomials form a system of complete orthogonal polynomials, which has numerous applications in science and engineering fields of study. The orthogonality can be expressed in the following way if \(P_{n} \left( x \right)\) denotes a Legendre polynomial of degree ‘\(n\)’

$$\mathop \int \limits_{ - 1}^{1} P_{m} \left( x \right)P_{n} \left( x \right)dx = 0\quad {\text{if}}\;n \ne m$$
(1)

Another distinguished property of the system is its definite parity, which is demonstrated by the relationship below

$$P_{n} \left( { - x} \right) = \left( { - 1} \right)^{n} P_{n} \left( x \right)$$
(2)

These properties make it convenient when Legendre polynomials are used in series expansion to approximate a function in the interval (− 1, 1). Also, the Legendre differential equation and the orthogonality property are independent of scale. The Legendre differential equation is given as

$$\left( {1 - x^{2} } \right)\frac{{d^{2} y}}{{dx^{2} }} - \frac{2xdy}{{dx}} + n\left( {n + 1} \right)y = 0$$
(3)

where \(n > 0\), \(\left| x \right| < 1\), or equivalently

$$\frac{d}{dx}\left[ {\left( {1 - x^{2} } \right)\frac{dy}{{dx}}} \right] + n\left( {n + 1} \right)y = 0$$
(4)

The solution of the above equations gives the Legendre functions of order ‘\(n\)’ with a general solution expressed as

$${\text{Y}} = {\text{A}}P_{n} \left( x \right) + BQ_{n} \left( x \right)$$
(5)

where \(P_{n} \left( x \right)\) and \(Q_{n} \left( x \right)\) are Legendre functions of the first and second kind of order ‘\(n\)’. For \(n\) = 1, 2, 3,…,\(P_{n} \left( x \right)\) functions are referred to as Legendre polynomials and are given by Rodrigue’s formula

$$P_{n} \left( x \right) = \frac{1}{{2^{n } n!}}\frac{{d^{n} }}{{dx^{n} }}\left( {x^{2} - 1} \right)^{n}$$
(6)

as well as in the recursive expression

$$\left( {n + 1} \right)P_{n + 1}^{^{\prime}} \left( x \right) = \left( {2n + 1} \right)xP_{n} \left( x \right) - nP_{n - 1}^{^{\prime}} \left( x \right)$$
(7)

3 The L-LSQ-FT and the L-IRLS-FT algorithms

As measured geophysical data always contain noise, the ability of a processing method to reduce or eliminate data noise is an advantageous feature. In order to add noise rejection capability to the Fourier transformation, we introduce the Legendre-Polynomials Least-Squares Fourier Transformation (L-LSQ-FT) and the Legendre-Polynomials Iteratively Reweighted Least-Squares Fourier Transformation (L-IRLS-FT).

3.1 Basic formulae for 1D case

In a one-dimensional case of physical applications, the Fourier transform (often abbreviated by FT) is well-defined as

$$U(\omega ) = \frac{1}{{\sqrt {2\pi } }}\int\limits_{ - \infty }^{\infty } {u(t)\,e^{ - j\omega t} dt} ,$$
(8)

where u represents a quantity related to a finite energy phenomenon and it is dependent on time denoted by \(t\), ω is the angular frequency and \(j\) symbolizes the imaginary unit. The frequency spectrum \(U(\omega )\) is a complex-valued continuous function representing the Fourier transform of a real-valued time function. Observably, a geophysical quantity measured in the time domain is adequately represented in the frequency domain by the Fourier transform. A backward transition from the frequency domain to the time domain is also possible and is provided by the inverse Fourier transform.

$$u(t) = \frac{1}{{\sqrt {2\pi } }}\int\limits_{ - \infty }^{\infty } {U(\omega )\,e^{j\omega t} d\omega }$$
(9)

Expressing the Fourier transform as an inverse problem required the frequency spectrum \(U(\omega )\) to be defined by a discrete parametric model. To realize this condition, we assumed that the parameter \(U(\omega )\) can be estimated with satisfactory accuracy if a finite series expansion is used as the approximation of the frequency spectrum

$$U(\omega ) \cong \sum\limits_{i = 1}^{M} {B_{i} \Psi_{i} (\omega )} ,$$
(10)

With the expression \(B_{i}\) being a complex-valued expansion coefficient and \(\Psi_{i}\) denoting a member of a properly chosen set of real-valued basis functions. The number of terms is limited to M. The theoretic value of a continuous signal in time domain at time moment tk can be given by the inverse Fourier transform as follows

$$u^{theor} (t_{k} ) = u_{k}^{theor} = \frac{1}{{\sqrt {2\pi } }}\int\limits_{ - \infty }^{\infty } {U(\omega )e^{{j\omega \,t_{k} }} d\omega } ,$$
(11)

Inserting the expression given in Eq. (9) one finds that

$$u_{k}^{theor} \cong \frac{1}{{\sqrt {2\pi } }}\int\limits_{a}^{b} {\left( {\sum\limits_{i = 1}^{M} {B_{i} \Psi_{i} (\omega )} } \right)} e^{{j\omega t_{k} }} d\omega = \sum\limits_{i = 1}^{M} {B_{i} \frac{1}{{\sqrt {2\pi } }}\int\limits_{a}^{b} {\Psi_{i} (\omega )} } e^{{j\omega t_{k} }} d\omega ,$$
(12)

where interval [a, b] is the domain of the selected set of basis functions. It is clear that the accuracy of this approximation is highly dependent on how the domain of basis functions correlates to the bandwidth of the continuous signal.

In practice, not a continuous but a discrete signal with finite number of samples is usually available for the estimation of frequency spectrum. Therefore, the frequency range is limited by the Nyquist-Shannon sampling theorem. So that this limited frequency range can be preserved, it must be mapped on the interval of basis functions by means of a linear scaling operation.

Introducing the notation ω’ for the scaled angular frequency, one can write that

$$G_{k,i} = \frac{1}{{\sqrt {2\pi } }}\int\limits_{a}^{b} {\Psi_{i} (\omega^{\prime } )e^{{j\omega^{\prime } t_{k} }} d\omega^{\prime } } ,$$
(13)

where \(G_{k,i}\) constitutes the elements of the Jacobian matrix of the size N-by-M. Here, N denotes the number of time domain samples. The model was parameterized by introducing the Legendre polynomials as basis functions to give

$$G_{k,n} = \frac{1}{{\sqrt {2\pi } }}\int\limits_{ - 1}^{1} {P_{n} (\omega^{\prime } )e^{{j\omega^{\prime } t_{k} }} d\omega^{\prime } }$$
(14)

In fact, Eq. (14) corresponds to the inverse Fourier transform of the Legendre polynomial of degree n. Thus, we can symbolize the general element of the Jacobian matrix in the following way

$$G_{k,n} = {\mathcal{F}}^{ - 1}_{k} \{ P_{n} (w^{\prime } )\}$$
(15)

The basic idea of introducing a new inversion-based Fourier Transformation method is to approximate the inverse Fourier transform of Pn(ω′) by means of the inverse discrete Fourier transform (IDFT) procedure:

$$G_{k,n} = IDFT_{k} \{ P_{n} (w^{\prime } )\} .$$
(16)

The values of the Pn(ω′) functions are accurate (noise-free), so the application of IDFT (or IFFT) is independent of the noise problem (of the data set). By using this procedure, the \(G_{k,n}\) elements of the Jacobian matrix can be generated numerically. The theoretical value of the discrete signal at time moment \(t_{k}\) is

$$u_{{}}^{theor} (t_{k} ) = u_{k}^{theor} = \sum\limits_{i = 1}^{M} {B_{n} G_{k,n} }$$
(17)

and the k-th element of the data deviation vector is written as

$$e_{k} = u_{k}^{meas} - u_{k}^{theor} = u_{k}^{meas} - \sum\limits_{n = 1}^{M} {B_{n}^{{}} \,G_{kn} } \,.$$
(18)

3.2 Basic formulae for 2D case

For a function u(x,y), the 2D Fourier transform can be calculated by the following integral

$$U(\omega_{x} ,\omega_{y} ) = \frac{1}{2\pi }\int\limits_{ - \infty }^{\infty } {} \int\limits_{ - \infty }^{\infty } {u\,(x,y)\,e^{{ - j\,(\omega_{x} x + \omega_{y} y)}} dx\,dy} .$$
(19)

The inverse transform is expressed by the formula below

$$u(x,y) = \frac{1}{2\pi }\int\limits_{ - \infty }^{\infty } {\,\int\limits_{ - \infty }^{\infty } {U(\omega_{x} } ,\,\omega_{y} )} \,e^{{j\,(\omega_{x} x + \omega_{y} y)}} d\omega_{x} d\omega_{y} ,$$
(20)

where U(ωx, ωy) is the 2D spatial-frequency spectrum, x and y represent the spatial coordinates and ωx, ωy denote the spatial angular frequencies. The expression below was used in discretizing the continuous Fourier spectrum

$$U(\omega_{x} ,\omega_{y} ) = \sum\limits_{n = 1}^{N} {\sum\limits_{m = 1}^{M} {B{}_{n,m}\Psi_{n} (\omega_{x} )\Psi_{m} } } (\omega_{y} ),$$
(21)

where \(\Psi_{n} (\omega_{x} )\) and \(\Psi_{m} (\omega_{y} )\) are frequency-dependent basis functions with the domain of interval [a, b]. The expansion coefficients are denoted by Bn,m, N and M are used for specifying the numbers of terms in the finite 2D series expansion. To develop the L-IRLS-FT algorithm in 2D, the same inversion procedures in the one-dimensional case was followed. The general form of the Jacobian matrix in the case of two-dimensional series expansion based inverse Fourier transform is given as

$$G_{nm}^{kl} = \frac{1}{{\sqrt {2\pi } }}\int\limits_{a}^{b} {\Psi_{n} (\omega_{x} )e^{{j\omega_{x} x_{k} }} d\omega_{x} } \cdot \frac{1}{{\sqrt {2\pi } }}\int\limits_{a}^{b} {\Psi_{m} (\omega_{y} )e^{{j\omega_{y} y_{l} }} d\omega_{y} } ,$$
(22)

where \(G_{nm}^{kl}\) is an element of the Jacobian matrix. Similarly to the 1D case, the linear scaling of independent variables cannot be avoided when the angular frequency ranges coming from the Nyquist-Shannon sampling theorem and the domain of basis functions do not match. Thus, the notations \(\omega_{x}^{^{\prime}}\) and \(\omega_{y}^{^{\prime}}\) must be introduced for the scaled spatial angular frequencies.

Parameterization of the model is achieved by introducing the Legendre polynomials as basis functions to give,

$$G_{nm}^{kl} = \frac{1}{{\sqrt {2\pi } }}\int\limits_{ - 1}^{1} {P_{n} (\omega^{\prime }_{x} )e^{{j\omega^{\prime }_{x} x_{k} }} d\omega^{\prime }_{x} } \cdot \frac{1}{{\sqrt {2\pi } }}\int\limits_{ - 1}^{1} {P_{m} (\omega^{\prime }_{y} )e^{{j\omega^{\prime }_{y} y_{l} }} d\omega^{\prime }_{y} }$$
(23)

or in a more formal notation

$$G_{n,m}^{kl} = {\mathcal{F}}_{k}^{ - 1} \left\{ {P_{n} \left( {\omega_{x}^{^{\prime}} } \right)} \right\} \cdot {\mathcal{F}}_{l}^{ - 1} \left\{ {P_{m} \left( {\omega_{y}^{^{\prime}} } \right)} \right\}$$
(24)

We employed the frequently used inverse discrete Fourier transform (IDFT) to calculate the inverse Fourier transform of Legendre polynomials,

$$G_{n,m}^{kl} = IDFT\left\{ {P_{n} \left( {\omega_{x}^{^{\prime}} } \right)} \right\} \cdot IDFT\left\{ {P_{m} \left( {\omega_{y}^{^{\prime}} } \right)} \right\}$$
(25)

By using this procedure the \(G_{n,m}^{kl}\) elements of the Jacobian matrix can be generated numerically. At this point, a new inversion method is to be defined. The theoretical value of the data at point (\(x_{k}\), \(y_{l}\)) is given as

$${\text{u}}^{theor} \left( {x_{k} ,y_{l} } \right) = u_{s}^{theor} = \mathop \sum \limits_{n = 1}^{N} \mathop \sum \limits_{m = 1}^{M} B_{n,m} G_{n,m}^{kl}$$

Which can be simplified further to

$$u_{s}^{theor} = \mathop \sum \limits_{i = 1}^{I} B_{i} G_{s,i}$$

(where \(i = n + \left( {m - 1} \right)N, s = k + \left( {l - 1} \right)k\))

The general element of the deviation vector can be given in the following form

$$e_{s} = u_{s}^{meas} - u_{s}^{theor} = u_{s}^{meas} - \sum\limits_{i = 1}^{I} {B_{i} \,G_{s,i} } \,.$$
(26)

3.3 The L-LSQ-FT and the robust L-IRLS-FT methods

Using Eqs. (18) and (26) we can calculate the L2-norm, the misfit function as

$$E_{2} = \sum\limits_{k = 1}^{N} {e_{k}^{2} = } \sum\limits_{k = 1}^{N} {\left( {u_{k}^{meas} - u_{k}^{theor} } \right)^{2} = } \sum\limits_{k = 1}^{N} {\left( {u_{k}^{meas} - \sum\limits_{n = 1}^{M} {B_{n} \,G_{kn} } } \right)^{2} } .$$

The minimization of this function gives the normal equation of the L-LSQ-FT (Gaussian Least Squares) method

$${\mathbf{G}}^{{\text{T}}} {\mathbf{G}}\overset{\lower0.5em\hbox{$\smash{\scriptscriptstyle\rightharpoonup}$}} {B} = {\mathbf{G}}^{{\text{T}}} \overset{\lower0.5em\hbox{$\smash{\scriptscriptstyle\rightharpoonup}$}} {u} ^{{meas}}$$

resulting in the solution

$$\overset{\lower0.5em\hbox{$\smash{\scriptscriptstyle\rightharpoonup}$}} {B} = \left( {{\mathbf{G}}^{{\text{T}}} {\mathbf{G}}} \right)^{{ - 1}} {\mathbf{G}}^{{\text{T}}} \overset{\lower0.5em\hbox{$\smash{\scriptscriptstyle\rightharpoonup}$}} {u} ^{{meas}}.$$

The Weighted Least Squares procedure (WLSQ) can be a robust inversion method with significant noise reduction capability and reduced outlier sensitivity. This is because each weighted data contributes to the solution according to its error margin, hence obtaining a good result (Dobroka et al. 2012). The deviation vector was minimized using Cauchy weights, and additionally reformed to Cauchy-Steiner weights. This was important because the Gaussian Least Squares Method (LSQ) is applicable only when data noise follows a regular distribution, and the algorithm minimizes the \(L_{2}\)-norm of the deviation vector between the theoretical and calculated data. Ideally, most geophysical measurements contain irregular noise with sparsely occurring outliers making the least-squares method (LSQ) less efficient for processing. The minimized weighted \(L_{2}\)-norm of the deviation vector is given as

$$E_{w} = \sum\limits_{k = 1}^{N} {w_{k} e_{k}^{2} }$$
(27)

Applying Steiner’s Most Frequent Value method (MFV), the scale parameter \(\varepsilon^{2}\) was determined from the data set in an internal iteration loop. The so-called Cauchy-Steiner weights are given by

$$w_{k} = \frac{{\varepsilon^{2} }}{{\varepsilon^{2} + e_{k}^{2} }},$$
(28)

where \(\varepsilon_{{^{{}} }}^{2}\) the Steiner’s scale factor called dihesion is determined iteratively.

That results in a non-linear inverse problem occasioned by an estimated misfit function which is non-quadratic due to the introduction of the Cauchy-Steiner weights (because \(e_{k}\) contains the unknown expansion coefficients). Practically, this is solvable by applying the method of the Iteratively Reweighted Least Squares (Scales et al. 1988) which linearizes the procedure. In a stepwise order of solution, a 0-th order solution \(\overset{\lower0.5em\hbox{$\smash{\scriptscriptstyle\rightharpoonup}$}} {B} ^{{(0)}}\) is derived by using the non-weighted LSQ method and the weights are recalculated as

$$w_{k}^{(0)} = \frac{{\varepsilon^{2} }}{{\varepsilon^{2} + (e_{k}^{(0)} )^{2} }}$$

with \(e_{k}^{(0)} = u_{k}^{measured} - u_{k}^{(0)}\), where \(u_{k}^{{(0)}} = \sum\limits_{{i = 1}}^{M} {B_{i}^{{(0)}} G_{{ki}} }\) and the expansion coefficients are given by the LSQ method. In the first iteration, the misfit function

$$E_{w}^{{(1)}} = \sum\limits_{{k = 1}}^{N} {w_{k}^{{(0)}} e_{k}^{{(1)^2}} }$$

is minimized resulting in the linear set of normal equations

$${\mathbf{G}}^{T} {\mathbf{W}}^{{(0)}} {\mathbf{G}}\vec{B}^{{(1)}} = {\mathbf{G}}^{T} {\mathbf{W}}^{{(0)}} \vec{u}^{{measured}}$$

The minimization of the new misfit function

$$E_{w}^{(2)} = \sum\limits_{k = 1}^{N} {w_{k}^{(1)} e_{k}^{(2)^2} }$$

gives \(\vec{B}^{(2)}\) which serves again for the calculation of \(w_{k}^{(2)} .\) This procedure is repeated giving the typical j-th iteration step

$${\mathbf{G}}^{T} {\mathbf{W}}^{(j - 1)} {\mathbf{G}}\,\vec{B}^{(j)} = {\mathbf{G}}^{T} {\mathbf{W}}^{(j - 1)} \vec{u}^{measured}$$
(29)

with the \({\mathbf{W}}^{(j - 1)}\) weighting matrix

$$W_{kk}^{(j - 1)} = w_{k}^{(j - 1)}$$
(30)

Every level of these iterations encloses an internal loop for the calculation of the Steiner’s scale parameter which is repetitive until an accurate stop condition is met.

3.3.1 Numerical testing in 1D

A time-domain signal (Fig. 1) was created to test the noise reduction capability of the newly developed methods, L-LSQ-FT, L-IRLS-FT and the traditional DFT in one dimension. The noiseless time function of the test data is given by the expression below

$$u(t\,) = \kappa {\kern 1pt} t^{\eta } e^{ - \lambda t} \sin (\mu t + \phi )$$
(31)

where the Greek letters represent the parameters of the signal. Specified fixed values for the signal parameters as follows: \(\kappa \approx 738.91\), \(\eta = {2}\), \(\lambda = {20}\), \(\mu = 40\pi\), \(\phi = \pi /4\).

Fig. 1
figure 1

Calculated noise-free waveform

Full size image

The noise-free waveform was sampled at regular intervals of 0.005 (sec) measurement points ranging over the time interval of [− 1, 1] and processed using the traditional DFT method to give both the real and imaginary parts of the noise-free Fourier spectrum (Fig. 2). The L-LSQ-FT spectrum (Fig. 3) was also calculated using Legendre polynomials of the (maximal) order of M = 300. For arithmetic purposes, the Fourier spectrum was calculated on data set converted between the interval [− 1, 1] in both x and y coordinates giving a suitable scale in the wavenumber domain. The output transformed spectra were similar for both the traditional DFT and the L-LSQ-FT. This demonstrates the applicability of L-LSQ-FT in processing noise-free data. Following the successful application of both methods to the noise-free signal, Gaussian and Cauchy noise were introduced into the noise-free signal (Fig. 1) for processing. In geophysical applications, the Gaussian noise distribution is occasionally encountered in data processing. The Gaussian noisy signal with 0 mean and 0.01 variance is given in Fig. 4.

Fig. 2
figure 2

Processed DFT spectrum of the noise-free data

Full size image
Fig. 3
figure 3

Processed L-LSQ-FT spectrum of the noise-free waveform

Full size image
Fig. 4
figure 4

The generated noisy signal with Gaussian noise

Full size image

On the other hand, random noise which does not follow a regular distribution across a survey area often occur. This type of noise is mostly introduced into survey data from external sources such as data acquisition or survey designs and equipment limitations. They are inherent in geophysical data and are not related to the subsurface body of interest. Random noise reduction is a valid step to increase the signal to noise ratio in geophysical applications with several methods developed over the years to achieve this purpose (Liu et al. 2006; Al-Dossary and Marfurt 2007; Liu et al. 2009). This includes the development of filters using various forms of transforms such as Wavelet Transform (Deighan and Watts 1997), S-Transform (Askari and Siahkooli 2008) and Fourier Transform (Dobroka et al. 2012). Failure to adequately suppress random noise affects the quality of processed data and interpretation.

Random noise following Cauchy distribution was added to the noise-free waveform to produce a noisy signal (Fig. 5) for processing. To demonstrate the noise reduction capability of the two methods, the Gaussian noisy signal (Fig. 4) was processed with the traditional DFT and the L-LSQ-FT methods. The resultant transformed spectra in the real and imaginary form are shown in Figs. 6 and 7 for the DFT and L-LSQ-FT respectively. The output signals show a considerable suppression of Gaussian noise by the L-LSQ-FT method compared to the traditional DFT method. Although a comparison between the real and imaginary spectrum as produced from the traditional DFT and the L-LSQ-FT does not show much improvement in output Fourier spectra in both methods, the L-LSQ-FT algorithm was able to reject a substantial amount of the introduced Gaussian noise.

Fig. 5
figure 5

The generated noisy signal with Cauchy noise

Full size image
Fig. 6
figure 6

Processed DFT spectrum of the Gaussian noisy signal (D = 1.03*10−2)

Full size image
Fig. 7
figure 7

Processed L-LSQ-FT spectrum of the Gaussian noisy signal (D = 8.2*10−3)

Full size image

For quantitative characterization of the results, we introduce the RMS distance between (a) and (b) data sets (for example noisy and noiseless) in the time domain (data distance)

$$d = \sqrt {\frac{1}{N}\sum\limits_{k = 1}^{N} {\left( {u^{(a)} (t_{k} ) - u^{(b)} (t_{k} )} \right)^{2} } } \,,$$
(32)

as well as the frequency domain (model or spectra distance)

$$D = \sqrt {\frac{1}{N}\sum\limits_{k = 1}^{N} {\left\{ {\left( {{\text{Re}} \left[ {U^{(a)} (f_{k} ) - U^{(b)} (f_{k} )} \right]} \right)^{2} + \left( {{\text{Im}} \left[ {U^{(a)} (f_{k} ) - U^{(b)} (f_{k} )} \right]} \right)^{2} } \right\}} } \;\,.$$
(33)

In the case of the Gaussian noise, the distance between the noisy and noiseless data sets, d = 0.1032. The model or spectra distance between the DFT spectrum of the noisy (contaminated with Gaussian noise) and the noiseless data sets gave D = 1.03*10−2. There is some improvement using the L-LSQ-FT method characterized by the spectra distance between the noiseless and the noisy spectra: D = 8.2*10−3 Similarly, the DFT gave a spectra distance D = 4.16*10−2 for spectrum produced from the noisy Cauchy signal whilst the L-LSQ-FT gave D = 2.43*10−2. From the analyses, a higher noise reduction capability was exhibited by the L-LSQ-FT method compared to the traditional DFT method. The results demonstrate the outlier and random noise sensitivity of the DFT and to some extent, the Least Squares Methods, hence there was the need to define a more robust method for outliers and random noise suppression. We, therefore, introduced the L-IRLS-FT Method. To test the noise reduction capability of the L-IRLS-FT, the Cauchy noisy signal (Fig. 5) was processed with the traditional 1D-DFT and L-IRLS-FT algorithms. The results are shown in Figs. 8 and 9 respectively. The comparison of output signals proves that the newly developed L-IRLS-FT algorithm was effective in reducing the Cauchy noise component of the noisy signal.

Fig. 8
figure 8

Processed DFT spectrum of the Cauchy noisy signal (D = 4.16*10−2)

Full size image
Fig. 9
figure 9

Processed L-IRLS-FT spectrum of the Cauchy noisy signal (D = 1.32*10−2

Full size image

To analytically examine the results, we applied the RMS distance between two data sets (for example noisy and noiseless) in the frequency domain as well as the model or spectra distance as indicated by Eqs. (32) and (33). The DFT gave a spectra distance D = 4.16*10−2 for spectrum produced from the noisy Cauchy signal whilst the L-IRLS-FT gave a spectra distance of D = 1.32*10−2. From the above analyses, the L-IRLS-FT method showed a high noise reduction capability when both regular and irregular noise was added to the noise-free waveform for processing. The outcome completely validates the outlier and random noise sensitivity of the traditional DFT method.

3.3.2 Numerical testing in 2D

The newly developed inversion-based Fourier transforms, 2D L-LSQ-FT and 2D L-IRLS–FT were tested on a noise-free data set. A rectangular test surface of size [− 1,1] units in both x and y directions was created (Fig. 10). An anomaly (u = 0.7) in the centre of size [− 0.2, 0.2] units in both directions was created with a noise-free regular background. The generated data were sampled at an interval of dx = dy = 0.02 units so the total quantity of data becomes N = 101*101. We then applied the 2D DFT and the 2D L-LSQ-FT algorithms to calculate the 2D Fourier spectrum of the noise-free discrete data set.

Fig. 10
figure 10

Generated noise-free 2D test surface

Full size image

Figures 11 and 12 show the absolute value (amplitude spectrum) produced by the traditional 2D DFT and the 2D L-LSQ-FT respectively. The 2D L-LSQ-FT spectrum was estimated from Legendre polynomials of the (maximal) order of M = 45. For scaling reasons, the data were initially changed between the interval [− 1, 1] in both x and y coordinates before the Fourier spectrums were calculated. Visual comparison of the two spectra shows approximately the same output, indicating the effectiveness of both methods in processing a noise-free dataset.

Fig. 11
figure 11

Processed 2D-DFT amplitude spectrum of the noise-free data with space frequencies (\(f_{x,} f_{y} )\) on the horizontal and vertical axes

Full size image
Fig. 12
figure 12

Processed 2D L-LSQ-FT amplitude spectrum of the noise-free data with space frequencies (\(f_{x,} f_{y} )\) on the horizontal and vertical axes

Full size image

After successfully applying both methods to the noise-free test surface, data noise following Gaussian and Cauchy distribution were introduced into the noise-free test surface for processing. The Gaussian noisy data is shown in Fig. 13 with a much rougher surface area. Random noise following Cauchy distribution was added to the test surface to produce outlier data with spikes (Fig. 14).

Fig. 13
figure 13

Generated 2D noisy signal with Gaussian noise

Full size image
Fig. 14
figure 14

Generated 2D noisy signal with Cauchy noise

Full size image

To demonstrate the noise reduction capability of the two methods, the Gaussian and Cauchy noisy data were processed with the traditional 2D DFT and the 2D L-IRLS-FT methods. To quantitatively describe the results, we introduce the RMS distance between (a) and (b) data in the space domain as

$$d_{RMS} = \sqrt {\frac{1}{N}\sum\limits_{i = 1}^{{N_{x} }} {\sum\limits_{j = 1}^{{N_{y} }} {\left[ {u^{(a)} (x_{i} ,y_{j} ) - u^{(b)} (x_{i} ,y_{j} )} \right] \, ^{2} } } }$$
(34)

And further in the space domain and the model distance as

$$D_{RMS} = \left[ \begin{gathered} \frac{1}{N}\sum\limits_{i = 1}^{{N_{x} }} {\sum\limits_{j = 1}^{{N_{y} }} {\left( {{\text{Re}} \, \left[ {U^{(a)} (\omega_{xi} ,\omega_{yj} )} \right] - {\text{Re}} \, \left[ {U^{(b)} (\omega_{xi} ,\omega_{yj} )} \right]} \right)^{2} } } \hfill \\ + \frac{1}{N}\sum\limits_{i = 1}^{{N_{x} }} {\sum\limits_{j = 1}^{{N_{y} }} {\left( {{\text{Im}} \, \left[ {U^{(a)} (\omega_{xi} ,\omega_{yj} )} \right] - {\text{Im}} \, \left[ {U^{(b)} (\omega_{xi} ,\omega_{yj} )} \right]} \right)}^{2} } \hfill \\ \end{gathered} \right]^{\frac{1}{2}} ,$$
(35)

(\(N_{x} ,N_{y} \,\) and \(N = N_{x} N_{y}\) are relevant numbers of data or space-frequency points in the 2D test area).

The output processed Fourier spectrum for the 2D DFT and the 2D L-IRLS-FT for the Cauchy noisy data are shown in Figs. 15 and 16 below. It appears from the output spectrums that the newly developed 2D L-IRLS-FT algorithm was effective in reducing the Cauchy noise component of the noisy data sets. The 2D L-IRLS-FT algorithm gave an almost noise-free spectrum demonstrating its effective noise reduction capability. For the processed Gaussian noisy dataset, the model or spectra distance between the DFT spectrum of the noisy (contaminated with Gaussian noise) and the noiseless DFT processed data sets was D = 1.700*10−3. There was sufficient improvement characterized by the spectra distance between the noiseless and the noisy (given by 2D L-IRLS-FT) spectra: D = 1.231*10−3. Likewise, the 2D DFT gave a spectra distance D = 3.413*10−3 for spectrum produced from the noisy Cauchy data (Fig. 15) whilst the 2D L-IRLS-FT gave an improved spectra distance of D = 6.759*10−4 (Fig. 16). Accordingly, the 2D L-IRLS-FT method showed a remarkable noise reduction capability when Cauchy noise was added to the test surface for processing. The results wholly validate the outlier and random noise sensitivity of the traditional 2D DFT method. Hence, we recommend the new method, the 2D L-IRLS-FT, when a robust and resistant solution is necessary for subduing the effects of random noise.

Fig. 15
figure 15

Processed 2D DFT spectrum of the Cauchy noisy Data with space frequencies (\(f_{x,} f_{y} )\) on the horizontal and vertical axes (D = 3.413*10−3)

Full size image
Fig. 16
figure 16

Processed 2D L-IRLS-FT spectrum of the Cauchy noisy Data with space frequencies (\(f_{x,} f_{y} )\) on the horizontal and vertical axes (D = 6.759*10−4)

Full size image

4 Application to reduction to the pole of magnetic data

The above tests prove the appreciable noise rejection capacity of the L-IRLS-FT method. This promising feature implies its application among others, in the field of reduction to the pole of magnetic data. As is well-known the shape and position of the maximum and minimum values of a magnetic anomaly are dependent on the inclination and declination angles. To eliminate this feature, the reduction to the magnetic pole is used. The result of this linear transformation is a magnetic anomaly of the same source assuming to be magnetized in the vertical direction (Kiss 1981). The transformation in the space-frequency \((f_{x} ,f_{y} )\,\) domain takes the form

$$T_{P} (f_{x} ,f_{y} ) = S_{P} (f_{x} ,f_{y} )\,T(f_{x} ,f_{y} )\,$$
(36)

where \(T(f_{x} ,f_{y} )\) is the Fourier transform of the total magnetic anomalies, \(S_{P} (f_{x} ,f_{y} )\,\) is the transfer function of the reduction to the pole operation

$$S_{P} (f_{x} ,f_{y} ) = (f_{x}^{2} + f_{y}^{2} )/(N(f_{x}^{2} + f_{y}^{2} )^{1/2} + j(Lf_{x} + Mf_{y} ))(n(f_{x}^{2} + f_{y}^{2} )^{1/2} + j(lf_{x} + mf_{y} )).$$

here L, M and l, m are direction cosines containing the declination and inclination angles (Kiss 1981). The total field reduced to pole is given by the inverse Fourier transformation

$$T_{P} (x,y) = {\mathscr {F}}_{k}^{{ - 1}} \left\{ {\left. {T_{P} (f_{x} ,f_{y} )} \right\}} \right..$$
(37)

The procedure contains the Fourier transform of the measured data, so in the case of noisy datasets containing outliers, it can be important to use FT algorithms with special noise rejection capacity. To demonstrate the efficiency of the L-IRLS-FT method (using Cauchy-Steiner weights) in the procedure of reduction to pole we defined a synthetic dataset calculated on a C-S shaped anomaly with uniform 1[A/m] magnetization in depth between (Al-Dossary and Marfurt 2007; Amundsen 1991) km. The magnetic field was calculated by using the method of Kuranatnam (1981). Figure 17a shows the noise-free magnetic field and Fig. 17b shows its pole reduced version calculated by the traditional DFT.

Fig. 17
figure 17

The noise-free magnetic field (a) and its reduced to the pole version using the traditional DFT (b)

Full size image

To simulate noisy dataset containing outliers 5% noise of Gaussian distribution was added to the noise-free data and a randomly selected 20% of the dataset was further contaminated by an additional 20% noise. The noisy dataset and its pole reduced version calculated by the traditional DFT is shown in Fig. 18a and b, respectively. As it can be seen the pole reduction results are very noisy.

Fig. 18
figure 18

The noisy magnetic field containing outliers (a) and its reduced to the pole version using the traditional DFT (b)

Full size image

Using the L-IRLS-FT method in the reduction to pole procedure of the same noisy dataset the result is shown in Fig. 19b. It can be seen that the new L-IRLS-FT algorithm using Cauchy-Steiner weights has appreciable outlier rejection also in the reduction to pole procedure.

Fig. 19
figure 19

The noisy magnetic field containing outliers (a) and its reduced to the pole version using the L-IRLS-FT method containing Cauchy-Steiner weights (b)

Full size image

5 Conclusion

We propose a new, advanced and vigorous inversion-based Fourier transformation algorithm, the Legendre-Polynomials Least-Squares Fourier Transformation (L-LSQ-FT) and Legendre-Polynomials Iteratively Reweighted Least-Squares Fourier Transformation (L-IRLS-FT) methods. The new algorithms use Legendre polynomials as basis functions for discretizing the resultant Fourier spectrum employing series expansion. In doing so, the expansion coefficients were estimated by resolving an over-determined inverse problem. The iteratively reweighted least squares method was introduced as a mean of robustifying the entire procedure. With this method, each data contributes to the solution according to its error margin, as individual data were weighted iteratively using the Cauchy-Steiner weights. The final algorithm was very tough and resistant to data outliers, especially sparsely distributed ones. This advantaged position was reached by comparing it with the traditional DFT algorithm in reduction to the pole of magnetic data. The L-IRLS-FT algorithm was found to be a very efficient, robust and resistant procedure with a higher noise rejection capability. Also, a test on synthetic data showed similar conclusions. Consequently, the newly developed L-IRLS-FT can be considered as a better alternative to the traditional DFT.