**Inversion-Based Fourier Transform as a New Tool for Noise Rejection** Inversion-Based Fourier Transform as a New Tool for Noise Rejection

Mihály Dobróka, Hajnalka Szegedi and Péter Vass Mihály Dobróka, Hajnalka Szegedi and

Additional information is available at the end of the chapter Additional information is available at the end of the chapter

http://dx.doi.org/10.5772/66338

### Abstract

Péter Vass

In this study, a new inversion method is presented for performing two-dimensional (2D) Fourier transform. The discretization of the continuous Fourier spectra is given by a series expansion with the scaled Hermite functions as square-integrable set of basis functions. The expansion coefficients are determined by solving an overdetermined inverse problem. In order to define a quick algorithm in calculating the Jacobian matrix of the problem, the special feature that the Hermite functions are eigenfunctions of the Fourier transformation is used. In the field of inverse problem theory, there are numerous procedures for noise rejection, so if the Fourier transformation is formulated as an inverse problem, these tools can be used to reduce the noise sensitivity. It was demonstrated in many case studies that the use of Cauchy-Steiner weights could increase the noise rejection capability of geophysical inversion methods. Following this idea, the two-dimensional Fourier transform is formulated as an iteratively reweighted least squares (IRLS) problem using Cauchy-Steiner weights. The new procedure is numerically tested using synthetic data.

Keywords: noise rejection in Fourier transformation, series expansion–based inversion, robust Fourier transformation, Hermite functions, reduction to pole

## 1. Introduction

In signal processing, the frequency spectrum of the time domain signals plays a very important role. In order to change over from the time domain to the frequency domain, the Fourier transform is applied most frequently. In the case of equidistantly sampled discrete time domain data sets, the so-called discrete Fourier transform (DFT) algorithm is used to

© The Author(s). Licensee InTech. This chapter is distributed under the terms of the Creative Commons © 2017 The Author(s). Licensee InTech. This chapter is distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/3.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

Attribution License (http://creativecommons.org/licenses/by/3.0), which permits unrestricted use, distribution, and eproduction in any medium, provided the original work is properly cited.

determine the discrete frequency spectrum. In the numerically very efficient Fast Fourier Transform algorithm (FFT), the spectrum is determined by solving a complete set of inhomogeneous linear algebraic set of equations.

The measured data set always contains noise, which is linearly projected into the frequency domain during Fourier transformation, so the traditional FT algorithms are sensitive to noise, most significantly to non-Gaussian one. On the other hand, it is well-known that in the framework of inverse problem theory there are a collection of methods with excellent noise rejection capability. For this reason, it was proposed to handle the 1D Fourier Transform as an overdetermined inverse problem [1].

In inverse problem theory, it is known that the simple least squares (LSQ) method gives optimal solution in case of Gaussian data noises while it is very sensitive for outliers. To reduce the effect of outlying data various (robust) inversion methods have been developed. The least absolute deviation (LAD) is one of the most frequently applied robust inversion method, which can be numerically realized by linear programing or by using the so-called iteratively reweighted least squares (IRLS) procedure [2]. In this case, the L1 norm of the deviation between the observed and predicted data is minimized. The IRLS procedure which iteratively recalculates the so-called Cauchy weights results in a very efficient robust inversion method [3]. In applying Cauchy inversion, the scale parameter of the Cauchy weights should be a priori known. This problem is solved in the framework of the most frequent value (MFV) method (developed by Steiner [4, 5]), where the scale parameter is derived from the data set. The weights given by the MFV method have been extensively used in various IRLS inversion problems. A successful application in joint inversion of seismic and geoelectric data was published by Dobróka et al. [6]. Szűcs et al. [7] reported a considerable improvement due to the use of Steiner's weights in the interpretation of borehole geophysical data. The Cauchy weights improved by Steiner's MFV method (the so-called Cauchy-Steiner weights) were successfully applied in robust tomography algorithms by Dobróka and Szegedi [8].

In previous papers by Szegedi and Dobróka [9], the 1D Fourier transformation was handled as robust inverse problem using IRLS algorithm with Cauchy-Steiner weights. It was shown that the noise sensitivity of the continuous Fourier transform (and its discrete variants DFT and FFT) was appropriately reduced by using robust inversion. Following a fruitful inversion strategy developed at the Geophysical Department of the University of Miskolc we used series expansion as a discretization tool. Series expansion–based inversion methods were successfully used in the interpretation of borehole geophysical data [10, 11] and also in processing induced polarization data [12]. In this study, we further develop the previously published inversion-based 1D Fourier transform algorithm by extending it to 2D cases.

## 2. Theoretical background for 1D algorithm

The Fourier transform and its inverse allow establishing a connection between the time and frequency domain. For the one-dimensional case the Fourier transform is defined as

Inversion-Based Fourier Transform as a New Tool for Noise Rejection http://dx.doi.org/10.5772/66338 5

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

where t denotes the time, ω is the angular frequency and j is the imaginary unit, UðωÞ is the Fourier transform of a real-valued time function uðtÞ. Thus, the Fourier transform provides the frequency domain representation of a phenomenon investigated by the measurement of some quantity in the time domain. By means of the inverse Fourier transform

$$
\mu(t) = \frac{1}{\sqrt{2\pi}} \int\_{-\infty}^{\infty} \mathcal{U}(\omega) \, e^{j\omega t} d\omega,\tag{2}
$$

we can return from the frequency domain to the time domain.

A next step in formulating the Fourier transform as an inverse problem is the discretization of the frequency spectrum UðωÞ. In order to satisfy this requirement, let us assume that UðωÞ is approximated with sufficient accuracy by using a finite series expansion

$$\mathcal{U}I(\omega) = \sum\_{i=1}^{M} B\_i \Psi\_i(\omega),\tag{3}$$

where the parameter Bi is a complex valued expansion coefficient and Ψ<sup>i</sup> is a member of an accordingly chosen set of real valued basis functions.

Using the terminology of (discrete) inverse problem theory, the theoretical values of time domain data (forward problem) can be given by the inverse Fourier transform

$$
\mu^{\rm thor}(t\_k) = \mu\_k^{\rm thor} = \frac{1}{\sqrt{2\pi}} \int\_{-\infty}^{\infty} \mathcal{U}(\omega) e^{j\omega \cdot t\_k} d\omega,\tag{4}
$$

where tk is the kth sampling time. Inserting the expression given in Eq. (3) one finds that

$$u\_k^{\text{ther}} \underset{\sqrt{2\pi}}{\mathbf{1}} \prod\_{i=1}^{\bullet} \left(\sum\_{i=1}^{M} B\_i \Psi\_i(\omega)\right) \mathbf{e}^{j\omega t\_k} d\omega = \sum\_{i=1}^{M} B\_i \frac{1}{\sqrt{2\pi}} \int\_{-\infty}^{\infty} \Psi\_i(\omega) \mathbf{e}^{j\omega t\_k} d\omega. \tag{5}$$

Let us introduce the notation

determine the discrete frequency spectrum. In the numerically very efficient Fast Fourier Transform algorithm (FFT), the spectrum is determined by solving a complete set of inhomo-

The measured data set always contains noise, which is linearly projected into the frequency domain during Fourier transformation, so the traditional FT algorithms are sensitive to noise, most significantly to non-Gaussian one. On the other hand, it is well-known that in the framework of inverse problem theory there are a collection of methods with excellent noise rejection capability. For this reason, it was proposed to handle the 1D Fourier Transform as an

In inverse problem theory, it is known that the simple least squares (LSQ) method gives optimal solution in case of Gaussian data noises while it is very sensitive for outliers. To reduce the effect of outlying data various (robust) inversion methods have been developed. The least absolute deviation (LAD) is one of the most frequently applied robust inversion method, which can be numerically realized by linear programing or by using the so-called iteratively reweighted least squares (IRLS) procedure [2]. In this case, the L1 norm of the deviation between the observed and predicted data is minimized. The IRLS procedure which iteratively recalculates the so-called Cauchy weights results in a very efficient robust inversion method [3]. In applying Cauchy inversion, the scale parameter of the Cauchy weights should be a priori known. This problem is solved in the framework of the most frequent value (MFV) method (developed by Steiner [4, 5]), where the scale parameter is derived from the data set. The weights given by the MFV method have been extensively used in various IRLS inversion problems. A successful application in joint inversion of seismic and geoelectric data was published by Dobróka et al. [6]. Szűcs et al. [7] reported a considerable improvement due to the use of Steiner's weights in the interpretation of borehole geophysical data. The Cauchy weights improved by Steiner's MFV method (the so-called Cauchy-Steiner weights) were successfully applied in robust tomography algorithms by Dobróka and

In previous papers by Szegedi and Dobróka [9], the 1D Fourier transformation was handled as robust inverse problem using IRLS algorithm with Cauchy-Steiner weights. It was shown that the noise sensitivity of the continuous Fourier transform (and its discrete variants DFT and FFT) was appropriately reduced by using robust inversion. Following a fruitful inversion strategy developed at the Geophysical Department of the University of Miskolc we used series expansion as a discretization tool. Series expansion–based inversion methods were successfully used in the interpretation of borehole geophysical data [10, 11] and also in processing induced polarization data [12]. In this study, we further develop the previously published

The Fourier transform and its inverse allow establishing a connection between the time and

frequency domain. For the one-dimensional case the Fourier transform is defined as

inversion-based 1D Fourier transform algorithm by extending it to 2D cases.

2. Theoretical background for 1D algorithm

geneous linear algebraic set of equations.

4 Fourier Transforms - High-tech Application and Current Trends

overdetermined inverse problem [1].

Szegedi [8].

$$\mathcal{G}\_{k,i} = \frac{1}{\sqrt{2\pi}} \int\_{-\infty}^{\infty} \Psi\_i(\omega) e^{j\omega t\_k} d\omega,\tag{6}$$

where Gk,<sup>i</sup> is an element of the so called Jacobian matrix of the size N-by-M (N is the number of time domain data and M is the number of unknown expansion coefficients). It is important for later considerations to note, that the Jacobian can be written as the inverse Fourier transform (in t ¼ tk) of the Ψ<sup>i</sup> basis function. The theoretical values take the linear form

$$
\mu\_k^{\text{theor}} = \sum\_{i=1}^{M} B\_i G\_{k,i}.\tag{7}
$$

The parameterization is always an important step in constructing an inversion algorithm. In Fourier transformation the frequency spectrum is defined over the interval (−∞,∞), so the set of basis functions should be defined in the same domain. In addition, the use of orthonormal functions for the series expansion is also proposed to the parameterization of the model. Because of these reasons we have chosen the set of scaled Hermite functions for discretization (their square-integrability ensures the existence of their Fourier transform).

If one tries to extend the concept of inversion-based Fourier transform for two-dimensional (2D) (or even multidimensional) case, a quick and simpler way of calculation can be advantageous. For this reason, consider the basic formulae of Hermite polynomials and Hermite functions.

The basic Hermite polynomials can be defined by the Rodriguez formula

$$h\_n^{(0)}(\\\omega) = (-1)^n e^{\omega^2} \left(\frac{d}{d\omega}\right)^n e^{-\omega^2}, \qquad \qquad n = 0, 1, 2, \ldots, \tag{8}$$

and also can be generated by the recursive formula

$$h\_{n+1}^{(0)}(\\\omega) = 2\omega \, h\_n^{(0)}(\\\omega) \text{--} 2 \, n \, h\_{n-1}^{(0)}(\\\omega),\tag{9}$$

where h ð0Þ <sup>0</sup> ðωÞ ¼ 1, h ð0Þ <sup>1</sup> ðωÞ ¼ 2ω. The Hermite polynomials fulfill the orthogonality condition

$$\int\_{-\infty}^{\infty} e^{-\omega^2} \cdot h\_n^{(0)}(\omega) \cdot h\_m^{(0)}(\omega) \, d\omega = 2^n n! \sqrt{\pi} \delta\_{nm}, \quad \delta\_{nm} = \begin{cases} 0, n \not\le m \\ 1, n = m \end{cases},\tag{10}$$

where δnm denotes the Kronecker symbol. Based on this formula, the basic Hermite functions can be defined as

$$H\_n^{(0)}(\\\omega) = \frac{e^{\frac{\omega^2}{2}} \cdot h\_n^{(0)}(\\\\\omega)}{\sqrt{\sqrt{\pi}} \, n! \, 2^n}, \qquad \qquad n = 0, 1, 2, \dots \tag{11}$$

Afterward the function Hð0<sup>Þ</sup> <sup>n</sup> ðωÞ is not only a complete orthogonal but an orthonormal system

$$\int\_{-\infty}^{\infty} H\_n^{(0)}(\omega) \cdot H\_m^{(0)}(\omega) \, d\omega = \delta\_{nm}, \quad \delta\_{nm} = \begin{cases} 0, n \not\le m \\ 1, n = m \end{cases} \tag{12}$$

There is an important special feature of Hermite functions, namely that they are the eigenfunctions of the Fourier transform [13]

$$\mathcal{F}\{H\_n^{(0)}(t)\} = (-j)^n H\_n^{(0)}(\omega),\tag{13}$$

and for the inverse Fourier transform, respectively

Inversion-Based Fourier Transform as a New Tool for Noise Rejection http://dx.doi.org/10.5772/66338 7

$$\mathcal{F}^{-1}\{H\_n^{(0)}(\omega)\}=(\mathfrak{j})^n H\_n^{(0)}(t). \tag{14}$$

As it was given in reference [14], the Hermite functions have to be modified by scaling because in geophysical applications the frequency covers wide ranges. The Rodriguez formula for modified Hermite polynomials takes the form

$$h\_n(\omega, \alpha) = (-1)^n e^{\alpha \cdot \alpha^2} \left(\frac{d}{d\omega}\right)^n e^{-\alpha \cdot \alpha^2},\tag{15}$$

and can be also generated by the recursive formula

utheor <sup>k</sup> ¼ ∑ M i¼1

(their square-integrability ensures the existence of their Fourier transform).

The basic Hermite polynomials can be defined by the Rodriguez formula

<sup>n</sup>þ<sup>1</sup>ðωÞ ¼ <sup>2</sup><sup>ω</sup> <sup>h</sup>ð0<sup>Þ</sup>

<sup>2</sup> � <sup>h</sup>ð0<sup>Þ</sup> <sup>n</sup> <sup>ð</sup>ω<sup>Þ</sup> ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ffiffiffi

<sup>F</sup>fHð0<sup>Þ</sup>

<sup>m</sup> <sup>ð</sup>ω<sup>Þ</sup> <sup>d</sup><sup>ω</sup> <sup>¼</sup> <sup>2</sup>nn! ffiffiffi

where δnm denotes the Kronecker symbol. Based on this formula, the basic Hermite functions

n e <sup>ω</sup><sup>2</sup> d dω � �<sup>n</sup> e −ω<sup>2</sup>

hð0<sup>Þ</sup>

6 Fourier Transforms - High-tech Application and Current Trends

ð0Þ

ð ∞

−∞ e −ω<sup>2</sup> � <sup>h</sup>ð0<sup>Þ</sup>

where h

ð0Þ

can be defined as

<sup>0</sup> ðωÞ ¼ 1, h

Afterward the function Hð0<sup>Þ</sup>

<sup>n</sup> ðωÞ¼ð−1Þ

h ð0Þ

<sup>n</sup> <sup>ð</sup>ωÞ � <sup>h</sup>ð0<sup>Þ</sup>

<sup>n</sup> <sup>ð</sup>ωÞ ¼ <sup>e</sup><sup>−</sup>ω<sup>2</sup>

<sup>n</sup> <sup>ð</sup>ωÞ � <sup>H</sup>ð0<sup>Þ</sup>

and also can be generated by the recursive formula

Hð0<sup>Þ</sup>

ð ∞

−∞ Hð0<sup>Þ</sup>

eigenfunctions of the Fourier transform [13]

and for the inverse Fourier transform, respectively

The parameterization is always an important step in constructing an inversion algorithm. In Fourier transformation the frequency spectrum is defined over the interval (−∞,∞), so the set of basis functions should be defined in the same domain. In addition, the use of orthonormal functions for the series expansion is also proposed to the parameterization of the model. Because of these reasons we have chosen the set of scaled Hermite functions for discretization

If one tries to extend the concept of inversion-based Fourier transform for two-dimensional (2D) (or even multidimensional) case, a quick and simpler way of calculation can be advantageous. For this reason, consider the basic formulae of Hermite polynomials and Hermite functions.

<sup>n</sup> <sup>ð</sup>ωÞ−<sup>2</sup> n h<sup>ð</sup>0<sup>Þ</sup>

<sup>1</sup> ðωÞ ¼ 2ω. The Hermite polynomials fulfill the orthogonality condition

<sup>π</sup> <sup>p</sup> <sup>δ</sup>nm, <sup>δ</sup>nm <sup>¼</sup> <sup>0</sup>, <sup>n</sup>≠<sup>m</sup>

<sup>π</sup> <sup>p</sup> <sup>n</sup>! <sup>2</sup> <sup>p</sup> <sup>n</sup> , <sup>n</sup> <sup>¼</sup> <sup>0</sup>, <sup>1</sup>, <sup>2</sup>, …: (11)

<sup>1</sup>, <sup>n</sup> <sup>¼</sup> <sup>m</sup> :

<sup>n</sup> ðωÞ, (13)

<sup>n</sup> ðωÞ is not only a complete orthogonal but an orthonormal system

�

<sup>m</sup> <sup>ð</sup>ω<sup>Þ</sup> <sup>d</sup><sup>ω</sup> <sup>¼</sup> <sup>δ</sup>nm, <sup>δ</sup>nm <sup>¼</sup> <sup>0</sup>, <sup>n</sup>≠<sup>m</sup>

n Hð0<sup>Þ</sup>

There is an important special feature of Hermite functions, namely that they are the

<sup>n</sup> ðtÞg ¼ ð−jÞ

�

BiGk,i: (7)

, n ¼ 0, 1, 2, :::, (8)

<sup>n</sup>−1ðωÞ, (9)

(10)

(12)

<sup>1</sup>, <sup>n</sup> <sup>¼</sup> <sup>m</sup> ,

$$h\_{n+1}(\omega, \alpha) = 2\omega\alpha \, h\_n(\omega, \alpha) \text{--} 2 \, n\alpha \, h\_{n-1}(\omega, \alpha), \tag{16}$$

where α is the scale factor and h0ðω, αÞ ¼ 1, h1ðω, αÞ ¼ 2αω [15]. The normalizing equation is

$$\int\_{-\alpha}^{\alpha} e^{-\alpha w^2} \cdot h\_n^{(0)}(\omega, \alpha) \cdot h\_m^{(0)}(\omega, \alpha) \, d\omega = \sqrt{\frac{\pi}{\alpha}} (2\alpha)^n n! \delta\_{nm}, \quad \delta\_{nm} = \begin{cases} 0, n \not\equiv m \\ 1, n = m \end{cases}.\tag{17}$$

Thus, the scaled Hermite functions can be defined as

$$H\_n(\omega, \alpha) = \frac{e^{-\frac{a}{\frac{\alpha}{2}}\omega^2} \cdot h\_n(\omega, \alpha)}{\sqrt{\sqrt{\frac{\pi}{\alpha}} n! \left(2\alpha\right)^n}}.\tag{18}$$

In this case the normalizing equation is

$$\int\_{-\boldsymbol{\sigma}}^{\boldsymbol{\sigma}} H\_n(\boldsymbol{\omega}, \boldsymbol{\alpha}) \cdot H\_m(\boldsymbol{\omega}, \boldsymbol{\alpha}) d\boldsymbol{\omega} = \delta\_{nm}, \quad \delta\_{nm} = \begin{cases} 0, n \not\models m \\ 1, n = m \end{cases} \tag{19}$$

Introducing the notation <sup>ω</sup><sup>0</sup> <sup>¼</sup> ffiffiffi <sup>α</sup> <sup>p</sup> <sup>ω</sup> the hnðω, <sup>α</sup><sup>Þ</sup> modified Hermite polynomials can be traced back to the hð0<sup>Þ</sup> <sup>n</sup> base polynomials. Substituting ω<sup>0</sup> into Eq. (15) we obtain

$$h\_n(\omega, a) = (\sqrt{a})^n (-1)^n e^{\omega^2} \left(\frac{d}{d\omega'}\right)^n e^{-\omega'^2} = (\sqrt{a})^n h\_n^{(0)}(\omega') = (\sqrt{a})^n h\_n^{(0)}(\sqrt{a}\omega). \tag{20}$$

Similarly, the modified Hermite function can also be traced back to the basic case (Hð0<sup>Þ</sup> <sup>n</sup> ). According to Eq. (18), we get the following formula

$$H\_n(\omega, a) = \frac{e^{\frac{\omega^2}{2}} (\sqrt{a})^n \, h\_n(\omega')}{\sqrt{\sqrt{\pi} \frac{1}{\sqrt{a}} \, n! \, 2^n a^n}} = \sqrt[4]{a} \frac{e^{\frac{\omega^2}{2}} \, h\_n(\omega')}{\sqrt{\sqrt{\pi}} \, n! \, 2^n} = \sqrt[4]{a} \, H\_n^{(0)}(\omega') = \sqrt[4]{a} \, H\_n^{(0)}(\sqrt{a}\omega). \tag{21}$$

Expanding the spectrum by means of the modified Hermite functions, in accordance with Eq. (6) the Jacobian matrix can be written as the inverse Fourier transform of the Hnðω, αÞ basis functions

### 8 Fourier Transforms - High-tech Application and Current Trends

$$\mathcal{G}\_{kn} = \frac{1}{\sqrt{2\pi}} \int\_{-\infty}^{\infty} H\_n(\omega, \alpha) \cdot e^{j\,\omega \cdot t} \, d\omega. \tag{22}$$

Using Eq. (21) one finds

$$G\_{\rm int} = \frac{1}{\sqrt{2\pi}} \int\_{-\infty}^{\infty} \sqrt[4]{\alpha} \, H\_{\rm n}^{(0)}(\omega') \cdot e^{i \,\omega \cdot t} \, d\omega,\tag{23}$$

or taking the notations ω t ¼ ω<sup>0</sup> t 0 , <sup>ω</sup><sup>0</sup> <sup>¼</sup> ffiffiffi <sup>α</sup> <sup>p</sup> <sup>ω</sup> and <sup>t</sup> <sup>0</sup> <sup>¼</sup> <sup>t</sup> ffiffi <sup>α</sup> <sup>p</sup> into account we have

$$\mathcal{G}\_{\rm kn} = \frac{1}{\sqrt[4]{a}} \frac{1}{\sqrt{2\pi}} \int\_{-\infty}^{\infty} H\_{\pi}^{(0)}(\omega') \cdot \left. \mathcal{e}^{i \ \omega' \ \ell} \right. \\ \left. \left. \mathcal{L} \omega' = \frac{1}{\sqrt[4]{a}} \right. \right| \left. \mathcal{F}^{-1} \{ H\_{\pi}^{(0)}(\omega') \} . \tag{24}$$

Using the properties of the base Hermite functions from Eq. (14) Eq. (24) can be rewritten in the following form

$$\mathcal{G}\_{\hbar n} = \frac{1}{\sqrt[4]{\alpha}} (j)^{\hbar} H\_{\hbar}^{(0)}(t') = \frac{1}{\sqrt[4]{\alpha}} (j)^{\hbar} H\_{\hbar}^{(0)} \left( \frac{t}{\sqrt{\alpha}} \right). \tag{25}$$

This is a very important result in further developing the inversion-based Fourier transform method because the Jacobian matrix can be produced quickly, as the procedure do not require integration. This is especially important in case of 2D (or higher dimensional) Fourier transform.

In accordance with Eq. (7) the theoretical data can be obtained as a linear expression of the expansion coefficients using the easily calculated elements of the Jacobian matrix. The general element of the deviation vector can be given in the following form

$$
\mathfrak{e}\_k = \mathfrak{u}\_k^{\text{measured}} \mathfrak{u}\_k^{\text{theor}}.\tag{26}
$$

In the framework of inverse problem theory, various methods are given for the minimization of appropriately chosen norm of the deviation vector resulting in an estimation of the expansion coefficients (Bestimated <sup>i</sup> ). After this, the real and imaginary part of the estimated spectrum can be calculated at any frequency as

$$\mathcal{U}^{\text{estimated}}(\omega) = \sum\_{i=1}^{M} \mathcal{B}\_i^{\text{estimated}} H\_i(\omega, \alpha). \tag{27}$$

### 3. Theoretical background for 2D algorithm

For the two-dimensional case the Fourier transform is defined as

Inversion-Based Fourier Transform as a New Tool for Noise Rejection http://dx.doi.org/10.5772/66338 9

$$\mathcal{U}(\omega\_x, \omega\_y) = \frac{1}{2\pi} \int\limits\_{-\infty}^{\infty} \int\limits\_{-\infty}^{\infty} \mu\left(\mathbf{x}, y\right) e^{\neg \left(\omega\_x x + \omega\_y y\right)} d\mathbf{x} \, dy,\tag{28}$$

where x, y denotes the spatial coordinates, ωx, ω<sup>y</sup> are the (angular) spatial frequencies and j is the imaginary unit. The frequency spectrum Uðωx, ωyÞ is the Fourier transform of a real valued function u ðx, yÞ and it is generally a complex valued continuous function. In two dimensions the forward problem giving the theoretical values of the space domain data can be defined by the two-dimensional inverse Fourier transform

$$u(\mathbf{x}, y) = \frac{1}{2\pi} \int\limits\_{-\infty}^{\infty} \int \mathcal{U}(\omega\_{\mathbf{x}}, \omega\_{y}) \, e^{j \left(\omega\_{\mathbf{x}}\mathbf{x} + \omega\_{y}y\right)} d\omega\_{\mathbf{x}} d\omega\_{y},\tag{29}$$

where Uðωx, ωyÞ denotes the 2D spatial frequency spectrum, which will be discretized using the scaled Hermite functions defined above

$$\mathcal{U}(\omega\_{\mathbf{x}}, \omega\_{\mathbf{y}}) = \sum\_{n=1}^{N} \sum\_{m=1}^{M} B\_{n, \, m} \, H\_n(\omega\_{\mathbf{x}}, \alpha) \, H\_m(\omega\_{\mathbf{y}}, \beta) \, \,, \tag{30}$$

where

Gkn <sup>¼</sup> <sup>1</sup>

Gkn <sup>¼</sup> <sup>1</sup>

0

ð ∞

−∞ Hð0<sup>Þ</sup> <sup>n</sup> ðω<sup>0</sup> Þ � e j ω0 t 0

ffiffiffiffiffiffi <sup>2</sup><sup>π</sup> <sup>p</sup>

, <sup>ω</sup><sup>0</sup> <sup>¼</sup> ffiffiffi

ð ∞

ffiffiffi <sup>α</sup> <sup>p</sup><sup>4</sup> <sup>H</sup>ð0<sup>Þ</sup> <sup>n</sup> ðω<sup>0</sup> Þ � e

<sup>α</sup> <sup>p</sup> <sup>ω</sup> and <sup>t</sup>

Using the properties of the base Hermite functions from Eq. (14) Eq. (24) can be rewritten in the

This is a very important result in further developing the inversion-based Fourier transform method because the Jacobian matrix can be produced quickly, as the procedure do not require integration. This is especially important in case of 2D (or higher dimensional) Fourier trans-

In accordance with Eq. (7) the theoretical data can be obtained as a linear expression of the expansion coefficients using the easily calculated elements of the Jacobian matrix. The general

In the framework of inverse problem theory, various methods are given for the minimization of appropriately chosen norm of the deviation vector resulting in an estimation of the expan-

> M i¼1

Bestimated

<sup>k</sup> −utheor

<sup>i</sup> ). After this, the real and imaginary part of the estimated spectrum

ek <sup>¼</sup> <sup>u</sup>measured

<sup>U</sup>estimatedðωÞ ¼ <sup>∑</sup>

3. Theoretical background for 2D algorithm

For the two-dimensional case the Fourier transform is defined as

<sup>0</sup> <sup>¼</sup> <sup>t</sup> ffiffi

<sup>d</sup>ω<sup>0</sup> <sup>¼</sup> <sup>1</sup>

ffiffiffi <sup>α</sup> <sup>p</sup><sup>4</sup> <sup>F</sup> �<sup>1</sup>

−∞

Using Eq. (21) one finds

following form

form.

sion coefficients (Bestimated

can be calculated at any frequency as

or taking the notations ω t ¼ ω<sup>0</sup> t

Gkn <sup>¼</sup> <sup>1</sup>

8 Fourier Transforms - High-tech Application and Current Trends

ffiffiffi α p4

1 ffiffiffiffiffiffi <sup>2</sup><sup>π</sup> <sup>p</sup>

Gkn <sup>¼</sup> <sup>1</sup>

ffiffiffi <sup>α</sup> <sup>p</sup><sup>4</sup> <sup>ð</sup>j<sup>Þ</sup> n Hð0<sup>Þ</sup> <sup>n</sup> ðt 0 Þ ¼ <sup>1</sup> ffiffiffi <sup>α</sup> <sup>p</sup><sup>4</sup> <sup>ð</sup>j<sup>Þ</sup> n Hð0<sup>Þ</sup> n

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

ffiffiffiffiffiffi <sup>2</sup><sup>π</sup> <sup>p</sup> ð ∞

−∞

Hnðω, αÞ � e

<sup>j</sup> <sup>ω</sup> <sup>t</sup> dω: (22)

<sup>j</sup> <sup>ω</sup> <sup>t</sup> dω, (23)

Þg: (24)

: (25)

<sup>k</sup> : (26)

<sup>i</sup> Hiðω, αÞ: (27)

<sup>α</sup> <sup>p</sup> into account we have

t ffiffiffi α p � �

<sup>f</sup>Hð0<sup>Þ</sup> <sup>n</sup> ðω<sup>0</sup>

$$H\_n(\omega\_x, \alpha) = \frac{e^{\frac{a \cdot a\_x^2}{2}} h\_n(\omega\_x, \alpha)}{\sqrt{\sqrt{\frac{\pi}{a}} n! \left(2\alpha\right)^n}}, \quad h\_n(\omega\_x, \alpha) = (-1)^n e^{a \cdot a\_x^2} \left(\frac{d}{d\omega\_x}\right)^n e^{-a \cdot a\_x^2},\tag{31}$$

$$H\_m(\omega\_y, \boldsymbol{\beta}) = \frac{e^{\frac{\beta^2 y^2}{2}} h\_m(\omega\_y, \boldsymbol{\beta})}{\sqrt{\sqrt{\frac{\pi}{\beta}}} m \, 1 \, (2\beta)^m} \,, \quad h\_m(\omega\_y, \boldsymbol{\beta}) = (-1)^m e^{\beta^2 \omega\_y^2} \left(\frac{d}{d\omega\_y}\right)^m e^{-\beta^2 \omega\_y^2} \, . \tag{32}$$

Using Eq. (29) the data calculated at the point ðxk, yl Þ

$$u(\mathbf{x}\_k, y\_l) = \frac{1}{2\pi} \int\_{-\infty}^{\infty} \int\_{-\infty}^{\infty} \sum\_{n=1}^{M} B\_{n, \, m} \, H\_n(\omega\_k, \alpha) \, H\_m(\omega\_\mathbf{y}, \beta) \, e^{j \left(\omega\_x \mathbf{x}\_k + \omega\_y \mathbf{y}\_l\right)} d\alpha\_\mathbf{x} d\alpha\_\mathbf{y} \,, \tag{33}$$

where k ¼ ð1, 2,…,KÞ, l ¼ ð1, 2, …, LÞdenote the sequence numbers of the measurement points along the x and y directions, respectively. By introducing the Jacobian matrix, we can write

$$\mu(\mathbf{x}\_k, \mathbf{y}\_l) = \sum\_{n=1}^{N} \sum\_{m=1}^{M} B\_{n,-m} \mathbf{G}\_{k,l}^{n,m},\tag{34}$$

where

$$\begin{split} G\_{k,l}^{n,m} &= \frac{1}{2\pi} \int\_{-\infty}^{\infty} \int\_{-\infty}^{\infty} H\_{\pi}(\omega\_{\mathbf{x}}, \alpha) \, H\_{m}(\omega\_{\mathbf{y}}, \beta) \, e^{j\left(\omega\_{\mathbf{x}}\mathbf{x}\_{k} + \omega\_{\mathbf{y}}\mathbf{y}\_{l}\right)} d\alpha\_{\mathbf{x}} d\alpha\_{\mathbf{y}} \\ &= \frac{1}{\sqrt{2\pi}} \int\_{-\infty}^{\infty} H\_{\pi}(\omega\_{\mathbf{x}}, \alpha) \, e^{j\left(\omega\_{\mathbf{x}}\mathbf{x}\_{k}\right)} d\alpha\_{\mathbf{x}} \, \frac{1}{\sqrt{2\pi}} \int\_{-\infty}^{\infty} H\_{m}(\omega\_{\mathbf{y}}, \beta) \, e^{j\omega\_{\mathbf{y}}\mathbf{y}\_{l}} d\alpha\_{\mathbf{y}} \,. \end{split} \tag{35}$$

Similar to Eq. (21)

$$H\_n(\omega\_x, \alpha) = \sqrt[4]{\alpha} \, H\_n^{(0)}(\sqrt{\alpha}\omega\_x) \,, \qquad H\_m(\omega\_{y'}\beta) = \sqrt[4]{\beta} \, H\_m^{(0)}(\sqrt{\beta}\omega\_y) \, \tag{36}$$

and the Jacobian takes the form

$$G\_{k,l}^{n,m} = \frac{\sqrt[4]{a\beta}}{\sqrt{2\pi}} \int\_{-\infty}^{\infty} H\_n^{(0)}(\sqrt{a}\omega\_x) \, e^{i\
u\_k \cdot \mathbf{x}\_k} d\omega\_x \, \frac{1}{\sqrt{2\pi}} \int\_{-\infty}^{\infty} H\_m^{(0)}(\sqrt{\beta}\omega\_y) \, e^{i\omega\_y y\_l} d\omega\_y \,. \tag{37}$$

Using the notations

$$
\boldsymbol{\omega}\_{\mathbf{x}} \cdot \mathbf{x}\_{k} = \boldsymbol{\omega}\_{\mathbf{x}}^{\prime} \cdot \mathbf{x}\_{k}^{\prime}, \ \boldsymbol{\omega}\_{\mathbf{x}}^{\prime} = \sqrt{\boldsymbol{\alpha}} \boldsymbol{\omega}\_{\mathbf{x}}, \ \mathbf{x}\_{k}^{\prime} = \frac{\mathbf{x}\_{k}}{\sqrt{\boldsymbol{\alpha}}}, \ \boldsymbol{\omega}\_{\mathbf{y}} \boldsymbol{y}\_{l} = \boldsymbol{\omega}\_{\mathbf{y}}^{\prime} \boldsymbol{\upbeta}\_{l}^{\prime}, \ \boldsymbol{\omega}\_{\mathbf{y}}^{\prime} = \sqrt{\boldsymbol{\beta}} \boldsymbol{\alpha}\_{\mathbf{y}}, \ \boldsymbol{y}\_{l}^{\prime} = \frac{\mathbf{y}\_{l}}{\sqrt{\boldsymbol{\beta}}}, \tag{38}
$$

we can write

$$\begin{split} G\_{k,l}^{n,m} &= \frac{1}{\sqrt[4]{a\beta\rho\sqrt{2\pi}}} \prod\_{-\infty}^{\infty} H\_{n}^{(0)}(\boldsymbol{\omega\_{x}}') \, e^{i\boldsymbol{\omega\_{x}}'\boldsymbol{\omega\_{x}'}} d\boldsymbol{\omega\_{x}}' \frac{1}{\sqrt{2\pi}} \prod\_{-\infty}^{\infty} H\_{m}^{(0)}(\boldsymbol{\omega\_{y}}') \, e^{i\boldsymbol{\omega\_{y}}'\boldsymbol{y\_{l}'}} d\boldsymbol{\omega\_{y}}' \\ &= \frac{1}{\sqrt[4]{a\beta}} \mathcal{F}^{-1}\{H\_{n}^{(0)}(\boldsymbol{\omega\_{x}}')\} \, \ \mathcal{F}^{-1}\{H\_{m}^{(0)}(\boldsymbol{\omega\_{y}}')\} \end{split} \tag{39}$$

and applying the well-known Eq. (14), the Jacobian matrix can be written in its final form (without integration)

$$G\_{k,l}^{n,m} = \frac{(\not j)^{n+m}}{\sqrt[4]{a\beta}} H\_n^{(0)}(\mathbf{x}\_k^{'}) \, H\_m^{(0)}(y\_l^{'}) \,. \tag{40}$$

The programming of the algorithm is quite simple after using the transformation of the indices i ¼ n þ ðm−1ÞN, s ¼ k þ ðl−1ÞK. With these notations, the total number of the unknown expansion coefficient is I ¼ N þ ðM−1ÞN ¼ NM and that of the measured data is S ¼ K þ ðL−1ÞK ¼ KL. The theoretical data can be calculated as

$$
\mu\_s^{\text{theor}} = \sum\_{i=1}^I B\_i \text{ G}\_{s,i}, \tag{41}
$$

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

Inversion-Based Fourier Transform as a New Tool for Noise Rejection http://dx.doi.org/10.5772/66338 11

$$e\_s = \mu\_s^{\text{measured}} - \sum\_{i=1}^{I} B\_i \ G\_{s,i} \tag{42}$$

with ði ¼ 1, …, I, s ¼ 1, …, SÞ. After this, the inverse problem can be formulated in a straightforward manner.

### 4. Inversion algorithms

G<sup>n</sup>,<sup>m</sup> <sup>k</sup>,<sup>l</sup> <sup>¼</sup> <sup>1</sup> 2π ð ∞

and the Jacobian takes the form

G<sup>n</sup>,<sup>m</sup> <sup>k</sup>,<sup>l</sup> ¼

> ′ xk ′ , ω<sup>x</sup>

G<sup>n</sup>,<sup>m</sup>

<sup>k</sup>,<sup>l</sup> <sup>¼</sup> <sup>1</sup> ffiffiffiffiffi

<sup>¼</sup> <sup>1</sup> ffiffiffiffiffi αβ <sup>p</sup><sup>4</sup> <sup>F</sup> �<sup>1</sup>

¼ KL. The theoretical data can be calculated as

Similar to Eq. (21)

Using the notations

ω<sup>x</sup> xk ¼ ω<sup>x</sup>

(without integration)

we can write

−∞

<sup>¼</sup> <sup>1</sup> ffiffiffiffiffiffi <sup>2</sup><sup>π</sup> <sup>p</sup> ð ∞

10 Fourier Transforms - High-tech Application and Current Trends

Hnðωx, <sup>α</sup>Þ ¼ ffiffiffi

ffiffiffiffiffi αβ p<sup>4</sup> ffiffiffiffiffiffi <sup>2</sup><sup>π</sup> <sup>p</sup> ð ∞

−∞

′ <sup>¼</sup> ffiffiffi

αβ p<sup>4</sup> ffiffiffiffiffiffi <sup>2</sup><sup>π</sup> <sup>p</sup> ð ∞

ð ∞

Hnðωx, αÞ Hmðωy, βÞ e

<sup>j</sup> <sup>ω</sup>xxkdω<sup>x</sup>

<sup>j</sup> <sup>ω</sup>xxkdω<sup>x</sup>

<sup>α</sup> <sup>p</sup> , <sup>ω</sup>yyl <sup>¼</sup> <sup>ω</sup><sup>y</sup>

<sup>f</sup>Hð0<sup>Þ</sup> <sup>m</sup> ðω<sup>y</sup> ′ Þg

<sup>n</sup> ðxk ′ <sup>Þ</sup> <sup>H</sup>ð0<sup>Þ</sup> <sup>m</sup> ðyl ′

and applying the well-known Eq. (14), the Jacobian matrix can be written in its final form

The programming of the algorithm is quite simple after using the transformation of the indices i ¼ n þ ðm−1ÞN, s ¼ k þ ðl−1ÞK. With these notations, the total number of the unknown expansion coefficient is I ¼ N þ ðM−1ÞN ¼ NM and that of the measured data is S ¼ K þ ðL−1ÞK

nþm ffiffiffiffiffi αβ <sup>p</sup><sup>4</sup> <sup>H</sup>ð0<sup>Þ</sup>

utheor <sup>s</sup> ¼ ∑ I i¼1

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

Hnðωx, αÞ e

<sup>n</sup> <sup>ð</sup> ffiffiffi

j ðωxxkþωyyl

−∞

1 ffiffiffiffiffiffi <sup>2</sup><sup>π</sup> <sup>p</sup> ð ∞

<sup>α</sup> <sup>p</sup> <sup>ω</sup>x<sup>Þ</sup> , Hmðωy, <sup>β</sup>Þ ¼ ffiffiffi

1 ffiffiffiffiffiffi <sup>2</sup><sup>π</sup> <sup>p</sup> ð ∞

−∞

′ yl ′ , ω<sup>y</sup>

−∞

Hð0<sup>Þ</sup> <sup>m</sup> ðω<sup>y</sup> ′ Þ e j ω<sup>y</sup> ′ yl ′ dω<sup>y</sup> ′

Hð0<sup>Þ</sup> <sup>m</sup> <sup>ð</sup> ffiffiffi <sup>β</sup> <sup>p</sup> <sup>ω</sup>y<sup>Þ</sup> <sup>e</sup>

Þ dωxdω<sup>y</sup>

Hmðωy, βÞ e

β p<sup>4</sup> Hð0<sup>Þ</sup>

′ <sup>¼</sup> ffiffiffi

<sup>β</sup> <sup>p</sup> <sup>ω</sup>y, yl

<sup>m</sup> <sup>ð</sup> ffiffiffi

<sup>β</sup> <sup>p</sup> <sup>ω</sup>yÞ, (36)

<sup>j</sup>ωyyldω<sup>y</sup> : (37)

′ <sup>¼</sup> yl ffiffiffi <sup>β</sup> <sup>p</sup> , (38)

Þ : (40)

Bi Gs,i, (41)

<sup>j</sup>ωyyldω<sup>y</sup> :

(35)

(39)

−∞

−∞

<sup>α</sup> <sup>p</sup><sup>4</sup> <sup>H</sup>ð0<sup>Þ</sup>

Hð0<sup>Þ</sup> <sup>n</sup> <sup>ð</sup> ffiffiffi <sup>α</sup> <sup>p</sup> <sup>ω</sup>x<sup>Þ</sup> <sup>e</sup>

<sup>α</sup> <sup>p</sup> <sup>ω</sup>x, xk

−∞

G<sup>n</sup>,<sup>m</sup> <sup>k</sup>,<sup>l</sup> <sup>¼</sup> <sup>ð</sup>j<sup>Þ</sup>

Hð0<sup>Þ</sup> <sup>n</sup> ðω<sup>x</sup> ′ Þ e j ω<sup>x</sup> ′ xk ′ dω<sup>x</sup> ′ 1 ffiffiffiffiffiffi <sup>2</sup><sup>π</sup> <sup>p</sup> ð ∞

<sup>f</sup>Hð0<sup>Þ</sup> <sup>n</sup> ðω<sup>x</sup> ′ Þg <sup>F</sup> �<sup>1</sup>

′ <sup>¼</sup> xk ffiffiffi If the measured data set contains Gaussian noise, the minimization of the L2 norm of the deviation vector is applied. This is the case of the least squares method when

$$E\_2 = \sum\_{k=1}^{N} e\_k^2 \tag{43}$$

is minimized resulting in the well-known set of the normal equations

$$\mathbf{G}^T \mathbf{G} \stackrel{\textstyle T}{B} = \mathbf{G}^T \stackrel{\text{measured}}{\textstyle}. \tag{44}$$

By solving these normal equations, we can give an estimate for the complex expansion coefficients, and both the real and imaginary parts of the LSQ estimated Fourier transform (LSQ-FT) can be calculated at any frequency by using

$$dL^{\text{estimated}}(\omega) = \sum\_{i=1}^{M} B\_i^{\text{estimated}} \Psi\_i(\omega). \tag{45}$$

As is well-known, the least squares method gives optimal results only when the data-noise follows Gaussian distribution. This distribution seldom occurs in practice so other norms of the deviation vector are introduced. In order to define a robust inversion algorithm, the minimization of the weighted norm

$$E\_w = \sum\_{k=1}^{N} w\_k e\_k^2. \tag{46}$$

with the so-called Cauchy weights

$$w\_k = \frac{\sigma^2}{\sigma^2 + e\_k^2}.\tag{47}$$

is suggested (here σ<sup>2</sup> is an accordingly chosen positive number). Using this norm for the solution of inverse problems provides reliable results even if the input data sets contain outliers [9].

There is a problem with inversion procedures involving Cauchy weights, namely the scale parameter should be a priori given. This difficulty can easily be solved by using Steiner weights [4]. In the framework of Steiner's most frequent value method, the scale parameter σ<sup>2</sup> is derived from data residuals in an internal iteration loop. In the (j + 1)th step of this procedure Steiner's scale factor ε<sup>2</sup> <sup>j</sup>þ<sup>1</sup> (called dihesion) can be calculated from ε<sup>2</sup> <sup>j</sup> as

$$\varepsilon\_{j+1}^2 = 3 \frac{\sum\_{k=1}^N \frac{e\_k^2}{\left(\varepsilon\_j^2 + e\_k^2\right)^2}}{\sum\_{s=1}^N \left(\frac{1}{\varepsilon\_j^2 + e\_s^2}\right)^2},\tag{48}$$

where the ε<sup>0</sup> starting value in the 0th step is given as

$$
\varepsilon\_0 \le \frac{\sqrt{3}}{2} \left( e\_{\text{max}} - e\_{\text{min}} \right). \tag{49}
$$

The stop criterion can be defined on an experimental basis (for example, a fixed number of iterations). After this the Cauchy weights are modified by using the (Steiner's) scale parameter (Cauchy-Steiner weights)

$$w\_k = \frac{\varepsilon^2}{\varepsilon^2 + e\_k^2}.\tag{50}$$

In the case of Cauchy-Steiner weights the misfit function given in Eq. (46) is nonquadratic (because ek contains the unknown expansion coefficients) and so the inverse problem is nonlinear which can be solved again by applying the method of the iteratively reweighted

least squares [2]. In the framework of this algorithm a 0th order solution B ⇀ð0Þ is derived by using the (nonweighted) LSQ method and the weights are calculated as

$$w\_k^{(0)} = \frac{\varepsilon^2}{\varepsilon^2 + (e\_k^{(0)})^2}.\tag{51}$$

with e ð0Þ <sup>k</sup> <sup>¼</sup> <sup>u</sup>measured <sup>k</sup> <sup>−</sup>u<sup>ð</sup>0<sup>Þ</sup> <sup>k</sup> , where <sup>u</sup><sup>ð</sup>0<sup>Þ</sup> <sup>k</sup> ¼ ∑ M i¼1 B<sup>ð</sup>0<sup>Þ</sup> <sup>i</sup> Gki and the expansion coefficients are given by the LSQ method. In the first iteration the misfit function

$$E\_w^{(1)} = \sum\_{k=1}^{N} w\_k^{(0)} e\_k^{(1)2} \tag{52}$$

is minimized resulting in the linear set of normal equations

$$\mathbf{G}^T \mathbf{W}^{(0)} \mathbf{G} \overrightarrow{B}^{(1)} = \mathbf{G}^T \mathbf{W}^{(0)} \overrightarrow{u}^{\text{measured}} \tag{53}$$

of the (linear) weighted least squares method where the Wð0<sup>Þ</sup> weighting matrix (independent of B !ð1<sup>Þ</sup> ) is of the diagonal form W<sup>ð</sup>0<sup>Þ</sup> kk <sup>¼</sup> <sup>w</sup><sup>ð</sup>0<sup>Þ</sup> <sup>k</sup> . Solving Eq. (53) one finds

Inversion-Based Fourier Transform as a New Tool for Noise Rejection http://dx.doi.org/10.5772/66338 13

$$\overrightarrow{\boldsymbol{B}}^{(1)} = (\mathbf{G}^T \mathbf{W}^{(0)} \mathbf{G})^{-1} \mathbf{G}^T \mathbf{W}^{(0)} \overrightarrow{\boldsymbol{u}}^{\text{measured}} \tag{54}$$

$$\mathfrak{u}\_k^{(1)} = \sum\_{i=1}^{M} B\_i^{(1)} \mathbf{G}\_{ki}, \quad \mathfrak{e}\_k^{(1)} = \mathfrak{u}\_k^{\text{measured}} \neg \mathfrak{u}\_k^{(1)}, \quad \mathfrak{w}\_k^{(1)} = \frac{\varepsilon^2}{\varepsilon^2 + (\mathfrak{e}\_k^{(1)})^2}. \tag{55}$$

The minimization of the new misfit function

[4]. In the framework of Steiner's most frequent value method, the scale parameter σ<sup>2</sup> is derived from data residuals in an internal iteration loop. In the (j + 1)th step of this procedure

<sup>j</sup> as

!<sup>2</sup> , (48)

<sup>2</sup> <sup>ð</sup>emax−eminÞ: (49)

: (50)

⇀ð0Þ

<sup>2</sup> : (51)

<sup>2</sup> (52)

<sup>i</sup> Gki and the expansion coefficients are given by the

is derived by

(53)

<sup>j</sup>þ<sup>1</sup> (called dihesion) can be calculated from ε<sup>2</sup>

∑ N k¼1

∑ N s¼1

ffiffiffi 3 p

The stop criterion can be defined on an experimental basis (for example, a fixed number of iterations). After this the Cauchy weights are modified by using the (Steiner's) scale parameter

> wk <sup>¼</sup> <sup>ε</sup><sup>2</sup> <sup>ε</sup><sup>2</sup> <sup>þ</sup> <sup>e</sup><sup>2</sup> k

In the case of Cauchy-Steiner weights the misfit function given in Eq. (46) is nonquadratic (because ek contains the unknown expansion coefficients) and so the inverse problem is nonlinear which can be solved again by applying the method of the iteratively reweighted

e2 k ðε2 <sup>j</sup> <sup>þ</sup> <sup>e</sup><sup>2</sup> k Þ 2

1 ε2 <sup>j</sup> þ e<sup>2</sup> s

ε2 <sup>j</sup>þ<sup>1</sup> ¼ 3

ε0≤

least squares [2]. In the framework of this algorithm a 0th order solution B

w<sup>ð</sup>0<sup>Þ</sup>

E<sup>ð</sup>1<sup>Þ</sup> <sup>w</sup> ¼ ∑ N k¼1 w<sup>ð</sup>0<sup>Þ</sup> <sup>k</sup> e ð1Þ k

GB !ð1<sup>Þ</sup>

<sup>¼</sup> <sup>G</sup><sup>T</sup>Wð0<sup>Þ</sup>

<sup>k</sup> . Solving Eq. (53) one finds

of the (linear) weighted least squares method where the Wð0<sup>Þ</sup> weighting matrix (independent of

u !measured

<sup>k</sup> <sup>¼</sup> <sup>ε</sup><sup>2</sup> ε<sup>2</sup> þ ðe ð0Þ <sup>k</sup> Þ

using the (nonweighted) LSQ method and the weights are calculated as

<sup>k</sup> ¼ ∑ M i¼1 B<sup>ð</sup>0<sup>Þ</sup>

where the ε<sup>0</sup> starting value in the 0th step is given as

12 Fourier Transforms - High-tech Application and Current Trends

Steiner's scale factor ε<sup>2</sup>

(Cauchy-Steiner weights)

<sup>k</sup> <sup>¼</sup> <sup>u</sup>measured

<sup>k</sup> <sup>−</sup>u<sup>ð</sup>0<sup>Þ</sup>

) is of the diagonal form W<sup>ð</sup>0<sup>Þ</sup>

<sup>k</sup> , where <sup>u</sup><sup>ð</sup>0<sup>Þ</sup>

is minimized resulting in the linear set of normal equations

G<sup>T</sup>Wð0<sup>Þ</sup>

kk <sup>¼</sup> <sup>w</sup><sup>ð</sup>0<sup>Þ</sup>

LSQ method. In the first iteration the misfit function

with e ð0Þ

B !ð1<sup>Þ</sup>

$$E\_w^{(2)} = \sum\_{k=1}^{N} w\_k^{(1)} e\_k^{(2)} \tag{56}$$

gives B !ð2<sup>Þ</sup> which serves again for the calculation of w<sup>ð</sup>2<sup>Þ</sup> <sup>k</sup> . This procedure is repeated giving the typical jth iteration step

$$\mathbf{G}^T \mathbf{W}^{(j-1)} \mathbf{G} \stackrel{\rightarrow}{B}^{(j)} = \mathbf{G}^T \mathbf{W}^{(j-1)} \stackrel{\text{measured}}{\boldsymbol{\mu}}^{\text{measured}} \tag{57}$$

with the Wðj−1<sup>Þ</sup> weighting matrix

$$\mathcal{W}\_{kk}^{(j-1)} = \mathfrak{w}\_k^{(j-1)}.\tag{58}$$

(Here we note that each step of these iterations contain an internal loop for the determination of the Steiner's scale parameter.) This iteration is repeated until a proper stop criterion is met.

### 5. Numerical investigations

In order to test our inversion-based Fourier transform we generated a 2D data set in a rectangular test area of the size [−1,1] units in both x and y directions (Figure 1). In the homogeneous background (with the theoretical model value u = 0), there is a rectangular anomaly (with u = 1.0) in the center of size [−0.2, 0.2] units in both directions. The sampling intervals were dx = dy = 0.04 units so the number of data is N = 51\*51. The 2D Fourier spectrum of the (noise-free) discrete data set was calculated by means of 2D DFT algorithm, Figure 2 shows its absolute value (amplitude spectrum).

To test the outlier sensitivity of the Fourier transformation algorithms, the noisy data set I was generated, in which random noise of Cauchy distribution (with 0 location and 0.02 scale parameters) were added to the noise-free data set shown in Figure 1. Data set I containing outliers is shown in Figure 3 and its DFT (amplitude) spectrum is shown in Figure 4. It can be seen that compared to Figure 2 the DFT spectrum is highly distorted proving a sufficient noise sensitivity to the traditional DFT.

For quantitative characterization of the results we introduce the RMS distance between two data sets (for example noisy and noiseless) as

Figure 1. The noise-free test surface.

Figure 2. The 2D amplitude spectrum of the noise-free data set calculated by DFT.

Figure 3. The noisy test surface.

Figure 1. The noise-free test surface.

14 Fourier Transforms - High-tech Application and Current Trends

Figure 2. The 2D amplitude spectrum of the noise-free data set calculated by DFT.

Figure 4. The 2D amplitude spectrum of the noisy data set calculated by DFT.

$$d = \sqrt{\frac{1}{N} \sum\_{i=1}^{N\_x} \sum\_{j=1}^{N\_y} [u^{\text{noiseless}}(\mathbf{x}\_i, y\_j) - u^{\text{noise}}(\mathbf{x}\_i, y\_j)]^2} \tag{59}$$

in the space domain (N, Nx, Nyare relevant numbers of data point in the 2D test area) and the model distance

$$D = \begin{bmatrix} \frac{1}{M} \sum\_{i=1}^{M\_x} \sum\_{j=1}^{M\_y} \left( \text{Re} [\boldsymbol{\mathcal{U}}^{\text{noise}}(\boldsymbol{\omega}\_{xi}, \boldsymbol{\omega}\_{yi})] - \text{Re} [\boldsymbol{\mathcal{U}}^{\text{noiseless}}(\boldsymbol{\omega}\_{xi}, \boldsymbol{\omega}\_{yi})] \right)^2 + \\\ \frac{M\_x}{M} \sum\_{i=1}^{M\_y} \sum\_{j=1}^{M\_y} \left( \text{Im} [\boldsymbol{\mathcal{U}}^{\text{noise}}(\boldsymbol{\omega}\_{xi}, \boldsymbol{\omega}\_{yi})] - \text{Im} [\boldsymbol{\mathcal{U}}^{\text{noiseless}}(\boldsymbol{\omega}\_{xi}, \boldsymbol{\omega}\_{yi})] \right)^2 \end{bmatrix} \tag{60}$$

in the spatial frequency domain (M, Mx, Myare relevant numbers of data points). The distance between the noisy and noiseless data sets is d = 0.0984. Using Eq. (60) we find the model distance between the DFT spectra of the noisy (contaminated with Cauchy noise) and the noiseless data sets: D = 0.0713.

If we apply our inversion based (IRLS-FT) method for the same noisy data set we get an estimated spectrum shown in Figure 5. Compared to the DFT spectrum (Figure 4) this figure represents sufficient improvement characterized by the model distance between the noiseless and the noisy (given by IRLS-FT) spectra: D = 0.00128.

Figure 5. The 2D amplitude spectrum of the noisy data set calculated by IRLS-FT.

Figure 6. The 2D inverse FT of the estimated spectrum.

It is well known that DFT and inverse DFT sequentially retrieve the noisy input data set exactly. In our inversion-based robust Fourier transform method we solve an overdetermined set of equations. In this case, it is important to see the space domain data set given by the inverse Fourier transform of the IRLS-FT spectrum. This is the so-called calculated data introduced previously in defining the IRLS-FT algorithm

$$
\mu\_s^{\text{theor}} = \sum\_{i=1}^I B\_i \ G\_{s,i}. \tag{61}
$$

The result is shown in Figure 6. Compared to the noisy data set, the new inversion-based Fourier transform method has appreciable noise rejection capability. This is characterized by the data distance between the noiseless data set and the space domain data calculated by the IRLS-FT method: d = 0.0140. It can be seen, that compared to the common DFT our inversionbased 2D Fourier transformation method has around 6–7 times lower noise sensitivity both in space domain and frequency domain.

### 6. Application

d ¼

model distance

D ¼

noiseless data sets: D = 0.0713.

1 <sup>M</sup> ∑ Mx i¼1 ∑ My j¼1 �

16 Fourier Transforms - High-tech Application and Current Trends

þ 1 <sup>M</sup> ∑ Mx i¼1 ∑ My j¼1 �

and the noisy (given by IRLS-FT) spectra: D = 0.00128.

Figure 5. The 2D amplitude spectrum of the noisy data set calculated by IRLS-FT.

1 <sup>N</sup> <sup>∑</sup> Nx i¼1 ∑ Ny j¼1

ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi

Re½Unoisyðωxi,ωyiÞ�−Re½Unoiselessðωxi,ωyiÞ�

Im½Unoisyðωxi,ωyiÞ�−Im½Unoiselessðωxi,ωyiÞ�

<sup>Þ</sup>−unoisyðxi, yj

vuut (59)

Þ�2

�2 þ

�2

1 2

(60)

<sup>½</sup>unoiselessðxi, yj

in the space domain (N, Nx, Nyare relevant numbers of data point in the 2D test area) and the

in the spatial frequency domain (M, Mx, Myare relevant numbers of data points). The distance between the noisy and noiseless data sets is d = 0.0984. Using Eq. (60) we find the model distance between the DFT spectra of the noisy (contaminated with Cauchy noise) and the

If we apply our inversion based (IRLS-FT) method for the same noisy data set we get an estimated spectrum shown in Figure 5. Compared to the DFT spectrum (Figure 4) this figure represents sufficient improvement characterized by the model distance between the noiseless

> The Fourier transformation is widely used in solving scientific or technical problems. Here we present a geophysical application in the field of processing geomagnetic data set. It is well known that the magnetic field has generally dipolar nature. It means that a magnetic body (i.e., wall fragments buried with soil in an archeo-geophysical measurement) usually produces

doubled anomaly (positive and negative) in the magnetic map depending on the geographical position of the measurement area. The only exceptions are the northern and the southern magnetic poles of the Earth and the magnetic equator. In order to simplify the interpretation of magnetic maps an elegant way was developed: the reduction to pole. This is a transformation resulting a magnetic data set that one would measure above the same magnetic body on the north (or southern) pole.

In order to apply our robust 2D IRLS-FT method a synthetic data set was generated. The measurement area was defined on the surface between (-100, 100) m in both of x and y direction. An anomaly of magnetization 100 nT (with D = 2.5° declination and I = 63° inclination) was assumed between the z-coordinates (20, 10) m. A rectangular measurement system was assumed with 5 m spacing in both directions (resulting in 1681 "measurement" data). The surface magnetic data calculated by means of the method of Kunaratnam [16] are shown in Figure 7. As it was mentioned, the interpretation of magnetic measurements is often supported by reducing the data to I = 90° pole. This can be done in the spatial frequency domain by applying the formula

$$R(\mu, \upsilon) = T(\mu, \upsilon) S(\mu, \upsilon), \tag{62}$$

where Tðu, vÞ is the 2D Fourier transform of the magnetic data set, Sðu, vÞ is the frequency domain operator of the pole reduction. The reduced data set in space domain can be found by inverse Fourier transformation of the Rðu, vÞ data set. This is shown in Figure 8 using noisefree magnetic data and the traditional DFT in 2D Fourier transformation.

Figure 7. The noise-free synthetic data set.

Figure 8. The data after reduction to pole (using DFT).

doubled anomaly (positive and negative) in the magnetic map depending on the geographical position of the measurement area. The only exceptions are the northern and the southern magnetic poles of the Earth and the magnetic equator. In order to simplify the interpretation of magnetic maps an elegant way was developed: the reduction to pole. This is a transformation resulting a magnetic data set that one would measure above the same magnetic body on

In order to apply our robust 2D IRLS-FT method a synthetic data set was generated. The measurement area was defined on the surface between (-100, 100) m in both of x and y direction. An anomaly of magnetization 100 nT (with D = 2.5° declination and I = 63° inclination) was assumed between the z-coordinates (20, 10) m. A rectangular measurement system was assumed with 5 m spacing in both directions (resulting in 1681 "measurement" data). The surface magnetic data calculated by means of the method of Kunaratnam [16] are shown in Figure 7. As it was mentioned, the interpretation of magnetic measurements is often supported by reducing the data to I = 90° pole. This can be done in the spatial frequency domain by

where Tðu, vÞ is the 2D Fourier transform of the magnetic data set, Sðu, vÞ is the frequency domain operator of the pole reduction. The reduced data set in space domain can be found by inverse Fourier transformation of the Rðu, vÞ data set. This is shown in Figure 8 using noise-

free magnetic data and the traditional DFT in 2D Fourier transformation.

Rðu, vÞ ¼ Tðu, vÞSðu, vÞ, (62)

the north (or southern) pole.

18 Fourier Transforms - High-tech Application and Current Trends

applying the formula

Figure 7. The noise-free synthetic data set.

Figure 9. The noisy synthetic data set.

Figure 10. The pole reduced data set (using DFT).

Figure 11. The pole reduced data set (using IRLS-FT).

In order to simulate noisy data set the magnetic data were contaminated with random noise following Cauchy distribution. The noisy data set and the result of pole reduction (using again the traditional DFT) is shown in Figures 9 and 10. It can be seen, that the pole reduction is highly distorted, which is caused by the low noise reduction capability of the 2D DFT algorithm proved in the previous chapter.

In contrary, the result of reduction to pole with the use of our new inversion-based 2D Fourier transformation algorithm is presented in Figure 11. In this case, we used Hermite function with 900 unknown expansion coefficients (considering the number of data, the inverse problem is sufficiently overdetermined). Figure 11 demonstrates high noise reduction capacity (compared to Figure 10, where the traditional 2D DFT was used for Fourier transformation). It can be seen that the pole-reduced data set is close to that shown in Figure 8 (noise-free data), and the limits of magnetization data are [0,250] in both cases. The result proves the successful applicability of our inversion-based 2D IRLS-FT algorithm.
