Efficient Application of the Voigt Functions in the Fourier Transform

1. Introduction

The forward and inverse Fourier transforms (FTs) can be defined in a symmetric form as [1,2]

$$\mathcal{F}\left\{f\left(t\right)\right\}\left(\nu \right)={\int }_{-\infty }^{\infty }f\left(t\right)e^{-2\pi i\nu t}dt=\widehat{f}\left(\nu \right)$$

(1)

and

$${\mathcal{F}}^{-1}\left\{\widehat{f}\left(\nu \right)\right\}\left(t\right)={\int }_{-\infty }^{\infty }\widehat{f}\left(\nu \right)e^{2\pi i\nu t}d\nu =f\left(t\right),$$

(2)

respectively. In this work we will consider only the forward FT (1) since due to symmetric form the approximations for the inverse FT (2) can be readily obtained from the forward FT by change of the variables.

There are relations between a function $f\left(t\right)$ and its even $f_{even}\left(t\right)$ and odd $f_{odd}\left(t\right)$ components. In particular, the even and odd functions can be readily generated by using the following relations

$$f_{even}\left(t\right)={\displaystyle \frac{f\left(t\right)+f(-t)}{2}}$$

(3)

and

$$f_{odd}\left(t\right)={\displaystyle \frac{f\left(t\right)-f(-t)}{2}}$$

(4)

such that

$$f\left(t\right)=f_{even}\left(t\right)+f_{odd}\left(t\right).$$

Therefore, due to linearity of the FT we can also state that

$$\widehat{f}\left(\nu \right)={\widehat{f}}_{even}\left(\nu \right)+{\widehat{f}}_{odd}\left(\nu \right).$$

(5)

Using Euler’s identity

$$e^{ix}=cos\left(x\right)+isin\left(x\right),$$

we can express equation (1) as

$$\mathcal{F}\left\{f\left(t\right)\right\}\left(\nu \right)={\int }_{-\infty }^{\infty }f\left(t\right)\left(cos\left(2\pi \nu t\right)-isin\left(2\pi \nu t\right)\right)dt$$

(6)

Since an even function satisfies the condition

$$f_{even}\left(t\right)=f_{even}(-t),$$

(7)

it follows that

$${\int }_{-\infty }^{\infty }{f}_{even}\left(t\right)cos\left(2\pi \nu t\right)dt=2{\int }_0^{\infty }{f}_{even}\left(t\right)cos\left(2\pi \nu t\right)dt$$

and

$${\int }_{-\infty }^{\infty }{f}_{even}\left(t\right)sin\left(2\pi \nu t\right)dt=0.$$

Since an odd function satisfies the condition

$$-f_{odd}\left(t\right)=f_{odd}(-t),$$

(8)

we have

$${\int }_{-\infty }^{\infty }{f}_{odd}\left(t\right)cos\left(2\pi \nu t\right)dt=0$$

and

$${\int }_{-\infty }^{\infty }{f}_{odd}\left(t\right)sin\left(2\pi \nu t\right)dt=2{\int }_0^{\infty }{f}_{odd}\left(t\right)sin\left(2\pi \nu t\right)dt.$$

Consequently, from these relations and equations (5), (6) we get the following identities

$$\mathcal{F}\left\{f_{even}\left(t\right)\right\}\left(\nu \right)=2{\int }_0^{\infty }{f}_{even}\left(t\right)cos\left(2\pi \nu t\right)dt={\widehat{f}}_{even}\left(\nu \right)$$

(9)

and

$$\mathcal{F}\left\{f_{odd}\left(t\right)\right\}\left(\nu \right)=-2i{\int }_0^{\infty }{f}_{odd}\left(t\right)sin\left(2\pi \nu t\right)dt={\widehat{f}}_{odd}\left(\nu \right).$$

(10)

Equations (9) and (10) are of the primary importance since they will be used in derivation of the new approximation of the FT.

In our recent publication we developed a new methodology providing a rational approximation of the FT [3]. Specifically, the FT of a function $f\left(t\right)$ can be approximated as a rational function consisting of low-order polynomials, $3\times 4$, in form

$$\begin{array}{cc}\hfill \widehat{f}\left(\nu \right)& \approx \sum _{m=1}^M\left({\displaystyle \frac{{\alpha }_m+{\beta }_m{\nu }^2}{{\nu }_m+{\lambda }_m{\nu }^2+{\nu }^4}}-i{\displaystyle \frac{{\gamma }_m\nu +{\theta }_m{\nu }^3}{{\kappa }_m+{\lambda }_m{\nu }^2+{\nu }^4}}\right)\hfill \\ \hfill & =\sum _{m=1}^M{\displaystyle \frac{{\alpha }_m-i{\gamma }_m\nu +{\beta }_m{\nu }^2-i{\theta }_m{\nu }^3}{{\kappa }_m+{\lambda }_m{\nu }^2+{\nu }^4}},\hfill \end{array}$$

(11)

where expansion coefficients are given by

$$\begin{array}{cc}\hfill {\alpha }_m& ={\displaystyle \frac{1}{8M{\pi }^4}}\sum _{n=-N}^N{f}_{even}\left(nh\right)e^{nh\sigma }\left({\mu }_m^2+{\sigma }^2\right)\left(\sigma cos\left(nh{\mu }_m\right)+{\mu }_{m}sin\left(nh{\mu }_m\right)\right),\hfill \\ \hfill {\beta }_m& ={\displaystyle \frac{1}{2M{\pi }^2}}\sum _{n=-N}^N{f}_{even}\left(nh\right)e^{nh\sigma }\left(\sigma cos\left(nh{\mu }_m\right)-{\mu }_{m}sin\left(nh{\mu }_m\right)\right),\hfill \end{array}$$

$$\begin{array}{cc}\hfill {\gamma }_m& ={\displaystyle \frac{1}{4M{\pi }^3}}\sum _{n=-N}^N{f}_{odd}\left(nh\right)e^{nh\sigma }\left({\sigma }^2-{\mu }_m^2\right)\left(cos\left(nh{\mu }_m\right)+2\sigma {\mu }_{m}sin\left(nh{\mu }_m\right)\right),\hfill \\ \hfill {\theta }_m& ={\displaystyle \frac{1}{M\pi }}\sum _{n=-N}^N{f}_{odd}\left(nh\right)e^{nh\sigma }cos\left(nh{\mu }_m\right),\hfill \\ \hfill {\kappa }_m& ={\displaystyle \frac{1}{16{\pi }^4}}{\left({\mu }_m^2+{\sigma }^2\right)}^2,\hfill \\ \hfill {\lambda }_m& ={\displaystyle \frac{1}{2{\pi }^2}}\left({\sigma }^2-{\mu }_m^2\right)\hfill \\ \hfill {\mu }_m& ={\displaystyle \frac{\pi (m-1/2)}{Mh}}.\hfill \end{array}$$

Some examples of rational approximation (11) of the FT are shown in our work [3] and can be visualized by running a provided Matlab code 1.

There are many methods for rational approximations such as Padé approximation [4,5], minimax approximation [6,7], Remez algorithm [8,9] and so on. However, to the best of our knowledge, a method of rational approximation that we developed in [3] is new and has been reported in scientific literature.

It should be noted that in our previous publication [10], we de facto applied a rational approximation of the FT. Therefore, a new method of rational approximation, described in our paper [3], is just a generalization of a sampling and integration technique that we proposed earlier in [10] to derive rapid and high-accuracy rational approximations of the complex error function.

The small order polynomials $3\times 4$ in quotients of rational approximation (11) may be advantageous for numerical analysis with any kind of computations including matrix manipulations, integrations and differentiations. However, the rational approximation (11) requires re-computation of the four expansion coefficients ${\alpha }_m$, ${\beta }_m$, ${\gamma }_m$, ${\theta }_m$ every time when the sampled function $f\left(t\right)$ changes. In this work, we propose a new method of rational approximation of the FT based on a sum of the real and imaginary parts of the complex error function [11,12,13]. Such a representation of rational approximation of the FT does not require re-computation when shape of the sampled function $f\left(t\right)$ changes.

2. Preliminaries

The complex error function, also commonly known as the Faddeeva function, can be defined as [13,14]

$$w\left(z\right)=e^{-z^2}\left(1+{\displaystyle \frac{2i}{\sqrt{\pi }}}{\int }_0^z{e}^{t^2}dt\right),$$

(12)

where $z=x+iy$. Comparing equation (12) with definition of the error function [11]

$$\mathrm{erf}\left(z\right)={\displaystyle \frac{2}{\sqrt{\pi }}}{\int }_0^z{e}^{-t^2}dt$$

one can see that the complex error function can be expressed as [11,15]

$$w\left(z\right)=e^{-z^2}\left[1-\mathrm{erf}(-iz)\right].$$

(13)

Therefore, the complex error function $w\left(z\right)$ can be considered as a reformulation of the error function $\mathrm{erf}\left(z\right)$.

Complex error function satisfies the following relation

$$w(-z)=2e^{-z^2}-w\left(z\right).$$

(14)

It can be shown that the complex error function is a solution of the following differential equation

$$w^{\prime }\left(z\right)+2zw\left(z\right)={\displaystyle \frac{2i}{\sqrt{\pi }}},$$

with initial condition

$$w\left(0\right)=1.$$

The complex error function $w\left(z\right)$ is closely related to the complex probability function [12,13]

$$W\left(z\right)={\displaystyle \frac{i}{\pi }}{\int }_{-\infty }^{\infty }{\displaystyle \frac{e^{-t^2}}{z-t}}dt.$$

(15)

There is a direct relationship between these two functions. In particular, both functions are equal to each other on the upper half of the complex plane [12,13]

$$w\left(z\right)=W\left(z\right),\mathrm{Im}\left[z\right]>0.$$

Separating the real and imaginary parts of the complex probability function (15) as

$$W\left(z\right)=K(x,y)+iL(x,y),$$

results in

$$K(x,y)={\displaystyle \frac{y}{\pi }}{\int }_{-\infty }^{\infty }{\displaystyle \frac{e^{-t^2}}{y^2+{(x-t)}^2}}dt$$

(16)

and

$$L(x,y)={\displaystyle \frac{1}{\pi }}{\int }_{-\infty }^{\infty }{\displaystyle \frac{e^{-t^2}(x-t)}{y^2+{(x-t)}^2}}dt,$$

(17)

respectively.

The real part $K(x,y)$ of the complex probability function is known as the Voigt function that is widely used in Atmospheric Physics to describe absorption and emission of atmospheric molecules [16,17,18].

The imaginary part $L(x,y)$ is also used in various fields of Physics and Engineering [19,20]. It does not have a specific name. However, following Zaghloul and Ali [15], for the convenience we will also refer to the function $L(x,y)$ as the imaginary Voigt function.

It is not difficult to show that substituting the following identity [12]

$${\displaystyle \frac{y}{y^2+{(x-t)}^2}}={\int }_0^{\infty }{e}^{-yq}cos\left[(x-t)q\right]dq,y>0,$$

(18)

into equation (16), we obtain [12,21]

$$\begin{array}{cc}\hfill K(x,y)& ={\displaystyle \frac{1}{\pi }}{\int }_{-\infty }^{\infty }{e}^{-t^2}{\displaystyle \frac{y}{y^2+{(x-t)}^2}}dt\hfill \\ & ={\displaystyle \frac{1}{\pi }}{\int }_0^{\infty }{\int }_{-\infty }^{\infty }{e}^{-t^2}{e}^{-yq}cos\left[(x-t)q\right]dtdq\hfill \\ & ={\displaystyle \frac{1}{\pi }}{\int }_0^{\infty }{e}^{-yq}{\int }_{-\infty }^{\infty }{e}^{-t^2}cos\left[(x-t)q\right]dtdq\hfill \end{array}$$

and since

$${\int }_{-\infty }^{\infty }{e}^{-t^2}cos\left[(x-t)q\right]dt=\sqrt{\pi }{e}^{-q^2/4}cos\left(qx\right),$$

we can write

$$K(x,y)={\displaystyle \frac{1}{\sqrt{\pi }}}{\int }_0^{\infty }{e}^{-t^2/4-yt}cos\left(xt\right)dt,y>0.$$

(19)

Similarly, substituting the identity

$${\displaystyle \frac{x-t}{y^2+{(x-t)}^2}}={\int }_0^{\infty }{e}^{-yq}sin\left(q(x-t)\right)dq,y>0,$$

into equation (17) we get [21]

$$\begin{array}{cc}\hfill L(x,y)& ={\displaystyle \frac{1}{\pi }}{\int }_{-\infty }^{\infty }{e}^{-t^2}{\displaystyle \frac{x-t}{y^2+{(x-t)}^2}}dt\hfill \\ & ={\displaystyle \frac{1}{\pi }}{\int }_0^{\infty }{\int }_{-\infty }^{\infty }{e}^{-t^2}{e}^{-yq}sin\left[(x-t)q\right]dtdq\hfill \\ & ={\displaystyle \frac{1}{\pi }}{\int }_0^{\infty }{e}^{-yq}{\int }_{-\infty }^{\infty }{e}^{-t^2}sin\left[(x-t)q\right]dtdq\hfill \end{array}$$

and since

$${\int }_{-\infty }^{\infty }{e}^{-t^2}sin\left[(x-t)q\right]dt=\sqrt{\pi }{e}^{-q^2/4}sin\left(qx\right),$$

we have

$$L(x,y)={\displaystyle \frac{1}{\sqrt{\pi }}}{\int }_0^{\infty }{e}^{-t^2/4-yt}sin\left(xt\right)dt,y>0.$$

(20)

Sum of the equations above in terms of the real and imaginary parts yields [21]

$$\begin{array}{cc}\hfill K(x,y)+iL(x,y)& ={\displaystyle \frac{1}{\sqrt{\pi }}}{\int }_0^{\infty }{e}^{-t^2/4-yt}\left[cos\left(xt\right)+isin\left(xt\right)\right]dt\hfill \\ & ={\displaystyle \frac{1}{\sqrt{\pi }}}{\int }_0^{\infty }{e}^{-t^2/4-yt}{e}^{ixt}dt\hfill \\ & ={\displaystyle \frac{1}{\sqrt{\pi }}}{\int }_0^{\infty }{e}^{-t^2/4}{e}^{(-y-ix)t}dt\hfill \\ & =e^{{(y-ix)}^2}\left[1-\mathrm{erf}(y-ix)\right].\hfill \end{array}$$

We can see now that this equation is consistent with equation (13) of the complex error function $w\left(z\right)$.

3. Results and Discussion

3.1. Methodology

Consider the rectangular function (solitary rectangular function) that can be defined as

$$f_r\left(t\right)=\left\{\begin{array}{cccc}& 1,\hfill & \hfill & -1/2<t<1/2,\hfill \\ & 1/2,\hfill & \hfill & \left|t\right|=1,\hfill \\ & 0,\hfill & \hfill & \mathrm{otherwise}.\hfill \end{array}\right.$$

This function can be expressed by the following limit

$$f_r\left(t\right)=\underset{k\to \infty }{lim}{\displaystyle \frac{1}{{\left(2t\right)}^{2k}+1}}.$$

(21)

Consequently, we can approximate the rectangular function by taking sufficiently large value of the parameter k.

Figure 1 shows the rectangular function and its approximation at $k=35$ by dashed red and light blue colors, respectively. As we can see, at $k=35$ equation (21) approximates the rectangular function reasonably well. Therefore, we can use the following approximation

$$f_r\left(t\right)\approx {\displaystyle \frac{1}{{\left(2t\right)}^{70}+1}}$$

(22)

and apply it to perform numerically the FT.

Figure 1. Rectangular and sawtooth functions and their approximations. Rectangular and sawtooth functions are shown by dashed red and solid black lines. Even $1/\left({\left(2t\right)}^{70}+1\right)$ and odd $t/\left({\left(2t\right)}^{70}+1\right)$ functions are shown by blue and green curves.

Figure 1. Rectangular and sawtooth functions and their approximations. Rectangular and sawtooth functions are shown by dashed red and solid black lines. Even $1/\left({\left(2t\right)}^{70}+1\right)$ and odd $t/\left({\left(2t\right)}^{70}+1\right)$ functions are shown by blue and green curves.

Previously, we have used the following sampling function

$$s\left(t\right)={\displaystyle \frac{h}{c\sqrt{\pi }}}{e}^{-{\left({\displaystyle \frac{t}{c}}\right)}^2},$$

where h and c are small fitting parameters, for high-accuracy approximation of the complex error function [22]. Thus, applying this sampling function over the points $nh$ to approximation (22), the rectangular function can be approximated as

$$\begin{array}{cc}\hfill f_r\left(t\right)& \approx \sum _{n=-N}^{N}s(t-nh)f_r\left(nh\right)\hfill \\ \hfill & ={\displaystyle \frac{h}{c\sqrt{\pi }}}\sum _{n=-N}^N{e}^{-{\left({\displaystyle \frac{t-nh}{c}}\right)}^2}{f}_r\left(nh\right)\hfill \\ \hfill & ={\displaystyle \frac{h}{c\sqrt{\pi }}}\sum _{n=-N}^N{e}^{\left(2{\displaystyle \frac{nh}{c^2}}\right)t}{e}^{-{\left({\displaystyle \frac{t}{c}}\right)}^2}{e}^{-{\left({\displaystyle \frac{nh}{c}}\right)}^2}{f}_r\left(nh\right)\hfill \end{array}$$

(23)

Figure 2 shows the approximations of the rectangular function at different fitting parameters h and c. Specifically, the green curve corresponds to $h=0.065$, $c=0.035$ and $N=25$, the red curve corresponds to $h=0.05$, $c=0.03$ and $N=25$ while the light blue curve corresponds to $h=0.02$, $c=0.025$ and $N=25$. The rectangular function is also shown by black dashed curve.

Since the rectangular function $f_r\left(t\right)$ is even, we can use equation (9) for the FT. Thus, the substitution of approximation (23) into equation (9) leads to

$$\begin{array}{cc}\hfill {\widehat{f}}_r\left(\nu \right)& \approx 2{\displaystyle \frac{h}{c\sqrt{\pi }}}{\int }_{-\infty }^{\infty }\sum _{n=-N}^N{e}^{-{\left({\displaystyle \frac{t}{c}}\right)}^2+\left(2{\displaystyle \frac{nh}{c^2}}\right)t}{e}^{-{\left({\displaystyle \frac{nh}{c}}\right)}^2}{f}_r\left(nh\right)cos(-2\pi \nu t)dt\hfill \\ & =2{\displaystyle \frac{h}{\sqrt{\pi }}}{\int }_{-\infty }^{\infty }\sum _{n=-N}^N{e}^{-t^2+{\displaystyle \frac{2nh}{c}}t}{e}^{-{\left({\displaystyle \frac{nh}{c}}\right)}^2}{f}_r\left(nh\right)cos\left(2\pi \nu ct\right)dt\hfill \\ & ={\displaystyle \frac{h}{\sqrt{\pi }}}\underset{-\infty }{\overset{\infty }{\int }}\sum _{n=-N}^N{e}^{-{\displaystyle \frac{t^2}{4}}+{\displaystyle \frac{nh}{c}}t}{e}^{-{\left({\displaystyle \frac{nh}{c}}\right)}^2}{f}_r\left(nh\right)cos\left(\pi \nu ct\right)dt.\hfill \end{array}$$

Comparing now approximation above with equation (19), one can see that the FT of the rectangular function (21) can be expressed in terms of the Voigt function

$${\widehat{f}}_r\left(\nu \right)=h\sum _{m=-N}^N{e}^{-{\left({\displaystyle \frac{nh}{c}}\right)}^2}{f}_r\left(nh\right)K\left(\pi \nu c,{\displaystyle \frac{nh}{c}}\right).$$

(24)

As a simplest example for an odd function, we can consider the sawtooth function (solitary sawtooth function)

$$f_s\left(t\right)=t{f}_r\left(t\right)=\left\{\begin{array}{cccc}& t,\hfill & \hfill & -1/2<t<1/2,\hfill \\ & 1/4,\hfill & \hfill & t=1,\hfill \\ & 0,\hfill & \hfill & \mathrm{otherwise}.\hfill \end{array}\right.$$

Since this function can be express though the following limit

$$f_s\left(t\right)=t\underset{k\to \infty }{lim}{\displaystyle \frac{1}{{\left(2t\right)}^{2k}+1}},$$

(25)

we can approximate it by multiplying t with equation (22) as follows

$$f_s\left(t\right)\approx t{\displaystyle \frac{1}{{\left(2t\right)}^{70}+1}}.$$

(26)

Figure 1 shows the sawtooth function (25) and its approximation (26) by black and light green curves, respectively.

Applying the sampling function over the points $nh$ to approximation (26), we obtain

$$\begin{array}{cc}\hfill f_s\left(t\right)& \approx \sum _{n=-N}^{N}s(t-nh)f_s\left(nh\right)\hfill \\ & ={\displaystyle \frac{h}{c\sqrt{\pi }}}\sum _{n=-N}^N{e}^{-{\left({\displaystyle \frac{t-nh}{c}}\right)}^2}{f}_s\left(nh\right)\hfill \\ & ={\displaystyle \frac{h}{c\sqrt{\pi }}}\sum _{n=-N}^N{e}^{\left(2{\displaystyle \frac{nh}{c^2}}\right)t}{e}^{-{\left({\displaystyle \frac{t}{c}}\right)}^2}{e}^{-{\left({\displaystyle \frac{nh}{c}}\right)}^2}{f}_s\left(nh\right).\hfill \end{array}$$

Consequently, the FT of the sawtooth function (25) can be approximated as

$$\begin{array}{cc}\hfill {\widehat{f}}_s\left(\nu \right)& \approx -2i{\displaystyle \frac{h}{c\sqrt{\pi }}}{\int }_{-\infty }^{\infty }\sum _{n=-N}^N{e}^{-{\left({\displaystyle \frac{t}{c}}\right)}^2+\left(2{\displaystyle \frac{nh}{c^2}}\right)t}{e}^{-{\left({\displaystyle \frac{nh}{c}}\right)}^2}{f}_s\left(nh\right)sin(-2\pi \nu t)dt\hfill \\ & =2i{\displaystyle \frac{h}{\sqrt{\pi }}}{\int }_{-\infty }^{\infty }\sum _{n=-N}^N{e}^{-t^2+{\displaystyle \frac{2nh}{c}}t}{e}^{-{\left({\displaystyle \frac{nh}{c}}\right)}^2}{f}_s\left(nh\right)sin\left(2\pi \nu ct\right)dt\hfill \\ & =i{\displaystyle \frac{h}{\sqrt{\pi }}}{\int }_{-\infty }^{\infty }\sum _{n=-N}^N{e}^{-{\displaystyle \frac{t^2}{4}}+{\displaystyle \frac{nh}{c}}t}{e}^{-{\left({\displaystyle \frac{nh}{c}}\right)}^2}{f}_s\left(nh\right)sin\left(\pi \nu ct\right)dt.\hfill \end{array}$$

Comparing this equation with the imaginary Voigt function (20), we can rearrange equation above as

$${\widehat{f}}_s\left(\nu \right)\approx ih\sum _{m=-N}^N{e}^{-{\left({\displaystyle \frac{nh}{c}}\right)}^2}{f}_s\left(nh\right)L\left(\pi \nu c,{\displaystyle \frac{nh}{c}}\right).$$

(27)

The applied methodology for the FTs of the rectangular and sawtooth functions (24) and (27) can be generalized to any even or odd functions as

$${\widehat{f}}_{even}\left(\nu \right)\approx h\sum _{m=-N}^N{e}^{-{\left({\displaystyle \frac{nh}{c}}\right)}^2}{f}_{even}\left(nh\right)K\left(\pi \nu c,{\displaystyle \frac{nh}{c}}\right)$$

(28)

and

$${\widehat{f}}_{odd}\left(\nu \right)\approx ih\sum _{m=-N}^N{e}^{-{\left({\displaystyle \frac{nh}{c}}\right)}^2}{f}_{odd}\left(nh\right)L\left(\pi \nu c,{\displaystyle \frac{nh}{c}}\right).$$

(29)

Due to symmetric properties of the even and odd functions (see equations (7) and (8)), the number of the summation terms in these approximations can be reduced by a factor of two.

Let us rearrange equation (28) in the following form

$$\begin{array}{cc}\hfill {\widehat{f}}_{even}\left(\nu \right)\approx & h\left[f_{even}\left(0\right)K(\pi \nu c,0)+\sum _{m=1}^N{e}^{-{\left({\displaystyle \frac{nh}{c}}\right)}^2}{f}_{even}\left(nh\right)K\left(\pi \nu c,{\displaystyle \frac{nh}{c}}\right)\right.\hfill \\ & \left.+\sum _{m=-N}^{-1}{e}^{-{\left({\displaystyle \frac{nh}{c}}\right)}^2}{f}_{even}\left(nh\right)K\left(\pi \nu c,{\displaystyle \frac{nh}{c}}\right)\right].\hfill \end{array}$$

Change of summation index as

$$\sum _{m=-N}^{-1}{e}^{-{\left({\displaystyle \frac{nh}{c}}\right)}^2}{f}_{even}\left(nh\right)K\left(\pi \nu c,{\displaystyle \frac{nh}{c}}\right)=\sum _{m=1}^N{e}^{-{\left({\displaystyle \frac{-nh}{c}}\right)}^2}{f}_{even}(-nh)K\left(\pi \nu c,{\displaystyle \frac{-nh}{c}}\right)$$

leads to

$$\begin{array}{cc}\hfill {\widehat{f}}_{even}\left(\nu \right)\approx & h\left[f_{even}\left(0\right)K(\pi \nu c,0)+\sum _{m=1}^N{e}^{-{\left({\displaystyle \frac{nh}{c}}\right)}^2}{f}_{even}\left(nh\right)K\left(\pi \nu c,{\displaystyle \frac{nh}{c}}\right)\right.\hfill \\ & \left.+\sum _{m=1}^N{e}^{-{\left({\displaystyle \frac{-nh}{c}}\right)}^2}{f}_{even}(-nh)K\left(\pi \nu c,{\displaystyle \frac{-nh}{c}}\right)\right].\hfill \end{array}$$

Since according to equation (7)

$$e^{-{\left({\displaystyle \frac{-nh}{c}}\right)}^2}=e^{-{\left({\displaystyle \frac{nh}{c}}\right)}^2},$$

$$f_{even}(-nh)=f_{even}\left(nh\right),$$

and since

$$K(x,0)=e^{-x^2}\Rightarrow K(\pi \nu c,0)=e^{-{\left(\pi \nu c\right)}^2},$$

we can write

$${\widehat{f}}_{even}\left(\nu \right)\approx h\left(f_{even}\left(0\right)e^{-{\left(\pi \nu c\right)}^2}+\sum _{m=1}^N{e}^{-{\left({\displaystyle \frac{nh}{c}}\right)}^2}{f}_{even}\left(nh\right)V_k\left(\pi \nu c,{\displaystyle \frac{nh}{c}}\right)\right),$$

(30)

where

$$V_k\left(\pi \nu c,{\displaystyle \frac{nh}{c}}\right)=K\left(\pi \nu c,{\displaystyle \frac{nh}{c}}\right)+K\left(\pi \nu c,{\displaystyle \frac{-nh}{c}}\right)$$

(31)

Equation (29) can be expanded as

$$\begin{array}{cc}\hfill {\widehat{f}}_{odd}\left(\nu \right)\approx & ih\left[f_{odd}\left(0\right)L(\pi \nu c,0)+\sum _{m=1}^N{e}^{-{\left({\displaystyle \frac{nh}{c}}\right)}^2}{f}_{odd}\left(nh\right)L\left(\pi \nu c,{\displaystyle \frac{nh}{c}}\right)\right.\hfill \\ \hfill & +\left.\sum _{m=-N}^{-1}{e}^{-{\left({\displaystyle \frac{nh}{c}}\right)}^2}{f}_{odd}\left(nh\right)L\left(\pi \nu c,{\displaystyle \frac{nh}{c}}\right)\right]\hfill \end{array}$$

(32)

Taking into consideration that

$$\sum _{m=-N}^{-1}{e}^{-{\left({\displaystyle \frac{nh}{c}}\right)}^2}{f}_{odd}\left(nh\right)L\left(\pi \nu c,{\displaystyle \frac{nh}{c}}\right)=\sum _{m=1}^N{e}^{-{\left({\displaystyle \frac{-nh}{c}}\right)}^2}{f}_{odd}(-nh)L\left(\pi \nu c,{\displaystyle \frac{-nh}{c}}\right)$$

such that according to (8)

$$f_{odd}(-nh)=-f_{odd}\left(nh\right)$$

and since

$$L(\pi \nu c,0)=0,$$

we can recast equation (32) as given by

$$\begin{array}{cc}\hfill {\widehat{f}}_{odd}\left(\nu \right)\approx & ih\left[\sum _{m=1}^N{e}^{-{\left({\displaystyle \frac{nh}{c}}\right)}^2}{f}_{odd}\left(nh\right)L\left(\pi \nu c,{\displaystyle \frac{nh}{c}}\right)\right.\hfill \\ & \left.+\sum _{m=1}^N{e}^{-{\left({\displaystyle \frac{-nh}{c}}\right)}^2}{f}_{odd}(-nh)L\left(\pi \nu c,{\displaystyle \frac{-nh}{c}}\right)\right]\hfill \end{array}$$

or

$$\begin{array}{cc}\hfill {\widehat{f}}_{odd}\left(\nu \right)\approx & ih[f_{odd}\left(0\right)L(\pi \nu c,0)\hfill \\ & +\sum _{m=1}^N{e}^{-{\left({\displaystyle \frac{nh}{c}}\right)}^2}{f}_{odd}\left(nh\right)\left(L\left(\pi \nu c,{\displaystyle \frac{nh}{c}}\right)-L\left(\pi \nu c,{\displaystyle \frac{-nh}{c}}\right)\right)]\hfill \end{array}$$

or

$${\widehat{f}}_{odd}\left(\nu \right)\approx ih\sum _{m=1}^N{e}^{-{\left({\displaystyle \frac{nh}{c}}\right)}^2}{f}_{odd}\left(nh\right)V_{\ell }\left(\pi \nu c,{\displaystyle \frac{nh}{c}}\right),$$

(33)

where

$$V_{\ell }\left(\pi \nu c,{\displaystyle \frac{nh}{c}}\right)=L\left(\pi \nu c,{\displaystyle \frac{nh}{c}}\right)-L\left(\pi \nu c,{\displaystyle \frac{-nh}{c}}\right).$$

(34)

At first glance the FT formulas (30) and (33) may appear computationally costly as they require user defined (external) function files. However, unlike equation (11), the expansion coefficients $V_k(\pi \nu c,nh/c)$ and $V_{\ell }(\pi \nu c,nh/c)$, shown by equations (31) and (34), respectively, do not require re-computations every time when we change the sampled function $f\left(t\right)$ to any other function. This gives a significant advantage. Since the values $V_k(\pi \nu c,nh/c)$ and $V_{\ell }(\pi \nu c,nh/c)$ remain always the same regardless the shape of the sampled function, these values can be precomputed in form of the look-up tables. Such implementation makes computation rapid as the required values $V_k(\pi \nu c,nh/c)$ and $V_{\ell }(\pi \nu c,nh/c)$ can be instantly picked up from the computer memory during computation of the FT. Furthermore, many algorithms, based on rational approximations, have been developed for rapid and high-accuracy computation of the Voigt/complex error function [23,24,25,26,27,28,29,30,31,32,33]. Therefore, the proposed technique can also be used as an alternative for rational approximation (11) of the FT.

3.2. Numerical Results

The FT of the rectangular function can be readily found by integration as

$${\widehat{f}}_r\left(\nu \right)={\int }_{-\infty }^{\infty }{f}_r\left(t\right)e^{-2\pi i\nu t}dt={\int }_{-1/2}^{1/2}{e}^{-2\pi i\nu t}dt=\mathrm{sinc}\left(\pi \nu \right),$$

(35)

where

$$\mathrm{sinc}\left(x\right)=\left\{\begin{array}{cccc}& {\displaystyle \frac{sin\left(x\right)}{x}},\hfill & \hfill & x\ne 0,\hfill \\ & 1,\hfill & \hfill & x=0,\hfill \end{array}\right.$$

is the sinc function [34,35]. Similarly, we can find the FT of the sawtooth function as follows

$$\begin{array}{cc}\hfill {\widehat{f}}_s\left(\nu \right)& ={\int }_{-\infty }^{\infty }{f}_s\left(t\right)e^{-2\pi i\nu t}dt={\int }_{-\infty }^{\infty }t{f}_r\left(t\right)e^{-2\pi i\nu t}dt\hfill \\ \hfill & ={\int }_{-1/2}^{1/2}t{e}^{-2\pi i\nu t}dt=i{\displaystyle \frac{\pi \nu cos\left(\pi \nu \right)-sin\left(\pi \nu \right)}{2{\pi }^2{\nu }^2}},\hfill \end{array}$$

(36)

We can use these analytical results for comparison with numerical FTs computed by using equations (30) and (33).

Figure 3 shows numerical FTs of rectangular and the sawtooth functions by light green and magenta curves, respectively. Equations (35) and (36) are also shown by dashed black curves.

Figure 4 illustrates the absolute differences ${\Delta }_r$ and ${\Delta }_s$ between equations (35), (36) and numerical FT approximations (30) and (33) at $h=0.02$, $c=0.025$ and $N=25$ by blue and red curves, respectively. As we can see, despite abrupt behavior of the rectangular and sawtooth functions $f_r\left(t\right)$ and $f_s\left(t\right)$ at $t_{1,2}=\pm 1/2$ and $t_{1,2}=\pm 1/4$, the FT approximations (24) and (27) can provide reasonable accuracies.

The accuracy of the FT can be significantly better for the well-behaved functions. As an example, consider a function

$$g\left(t\right)=e^{-{\left(6\pi t\right)}^2}-sin\left(32t\right)e^{-{\left(7\pi t\right)}^2}.$$

(37)

The first term of this function is even

$$g_{even}\left(t\right)={\displaystyle \frac{g\left(t\right)+g(-t)}{2}}=e^{-{\left(6\pi t\right)}^2}$$

(38)

while its second term is odd

$$g_{odd}\left(t\right)={\displaystyle \frac{g\left(t\right)-g(-t)}{2}}=-sin\left(32t\right)e^{-{\left(7\pi t\right)}^2}.$$

(39)

Figure 5 shows the functions $g\left(t\right)$, $g_{even}\left(t\right)$ and $g_{odd}\left(t\right)$ by blue, red and green colors, respectively.

The FTs of even and odd components (38) and (39) can be obtained analytically. In particular, substituting equations (38) and (39) into FT formulas (9) and (10) yields

$${\widehat{g}}_{even}\left(\nu \right)={\displaystyle \frac{e^{-{\left(\frac{\nu }{6}\right)}^2}}{6\sqrt{\pi }}}$$

(40)

and

$${\widehat{g}}_{odd}\left(\nu \right)=i{\displaystyle \frac{1}{14\sqrt{\pi }}}{e}^{-{\left({\displaystyle \frac{16+\pi \nu }{7\pi }}\right)}^2}\left(e^{{\displaystyle \frac{64\nu }{49\pi }}}-1\right).$$

(41)

Figure 6 depicts the absolute errors ${\Delta }_{even}$ and ${\Delta }_{odd}$ between equations (40), (41) and corresponding FT approximations (30) and (33) at $h=0.004$, $c=0.0045$ and $N=30$ by blue and red colors, respectively. As we can see, the absolute errors do not exceed $0.00035$ and $0.0005$ for the functions and ${\widehat{g}}_{even}\left(\nu \right)$ and ${\widehat{g}}_{odd}\left(\nu \right)$, respectively, over the wide interval of $\nu$.

3.3. Trigonometric Forms

As it has been mentioned above, the equations (28), (29) and their variations (30), (33) can be used to implement a rational approximation of the FT. We can also show how to represent these equations in trigonometric forms.

The real part of equation (14) gives

$$K(x,-y)=2e^{-x^2+y^2}cos\left(2xy\right)-K(x,y)$$

Therefore, equation (31) can be further simplified as

$$\begin{array}{cc}\hfill V_k\left(\pi \nu c,{\displaystyle \frac{nh}{c}}\right)& =K\left(\pi \nu c,{\displaystyle \frac{nh}{c}}\right)+2e^{-{\left(\pi \nu c\right)}^2+{\displaystyle \frac{nh}{c}}}cos\left(\pi \nu c,{\displaystyle \frac{nh}{c}}\right)-K\left(\pi \nu c,{\displaystyle \frac{nh}{c}}\right)\hfill \\ & =2e^{-{\left(\pi \nu c\right)}^2+{\displaystyle \frac{nh}{c}}}cos\left(\pi \nu c,{\displaystyle \frac{nh}{c}}\right).\hfill \end{array}$$

leading to

$${\widehat{f}}_{even}\left(\nu \right)\approx h{e}^{-{\left(\pi \nu c\right)}^2}\left(f_{even}\left(0\right)+2\sum _{m=1}^N{f}_{even}\left(nh\right)cos\left(2\pi \nu nh\right)\right)$$

(42)

in accordance with equation (30).

The imaginary part of equation (14) provides

$$L(x,-y)=2e^{-x^2+y^2}sin\left(2xy\right)+L(x,y)$$

Substituting this equation into approximation (34) yields

$$\begin{array}{cc}\hfill V_{\ell }\left(\pi \nu c,{\displaystyle \frac{nh}{c}}\right)& =L\left(\pi \nu c,{\displaystyle \frac{nh}{c}}\right)-2e^{-{\left(\pi \nu c\right)}^2+{\left({\displaystyle \frac{nh}{c}}\right)}^2}sin\left(2\pi \nu nh\right)-L\left(\pi \nu c,{\displaystyle \frac{nh}{c}}\right)\hfill \\ & =-2e^{-{\left(\pi \nu c\right)}^2+{\left({\displaystyle \frac{nh}{c}}\right)}^2}sin\left(2\pi \nu nh\right).\hfill \end{array}$$

This results in

$${\widehat{f}}_{odd}\left(\nu \right)\approx -2ih{e}^{-{\left(\pi \nu c\right)}^2}\sum _{m=1}^N{f}_{odd}\left(nh\right)\left(sin\left(2\pi \nu nh\right)\right)$$

(43)

according to equation (33).

As we can see, equations (42) and (43) represent trigonometric versions of the equations (30) and (33), respectively. Despite trigonometric representations, equations (42) and (43) do not have any advantage over equations (30) and (33) in computational speed since all required values $V_k(\pi \nu c,nh/c)$ and $V_{\ell }(\pi \nu c,nh/c)$ can be precomputed and saved in a computer memory in form of the look-up tables.

Combining equations (3), (4), (5), (42) and (43) together, we obtain

$$\begin{array}{cc}\hfill \widehat{f}\left(\nu \right)\approx & h{e}^{-{\left(\pi \nu c\right)}^2}[f\left(0\right)+\sum _{m=1}^N\left(\left(f\left(nh\right)+f(-nh)\right)cos\left(2\pi \nu nh\right)\right.\hfill \\ & \left.-i\left(f\left(nh\right)-f(-nh)\right)\left(sin\left(2\pi \nu nh\right)\right)\right)]\hfill \end{array}$$

or

$$\widehat{f}\left(\nu \right)\approx h{e}^{-{\left(\pi \nu c\right)}^2}\left(f\left(0\right)+\sum _{m=1}^N\left(f\left(nh\right)e^{-2\pi i\nu nh}+f(-nh)e^{2\pi i\nu nh}\right)\right).$$

(44)

We can see that equation (44) resembles the discrete Fourier transfer (DFT) [1,2]. However, unlike the DFT, equation (44) provides non-periodic output due to exponential multiplier $e^{-{\left(\pi \nu c\right)}^2}$ acting like the Hamming window that sometimes may be introduced to eliminate periodicity in time or frequency domains. Remarkably, this exponential multiplier in equation (44) was not introduced but appeared naturally in derivation.

Examples of the numerical FTs based on equations (30) and (33) can be validated by running the Mathematica codes, shown in the next section.

4. Mathematica Codes

4.1. Complex Error Function

To compute the complex error function $w\left(z\right)$, we can use, for example, the algorithm provided in our work [29]. To cover the entire complex plane with high accuracy of computation, we used only three approximations. We will use four cells to define each equation and to distribute them accordingly on the complex plane.

The first cell below is to define Equation (7) from [29].

Clear[a0, b0, c0, \[CapitalOmega], w1];

(* Fitting parameters *)

H = 0.25; \[Stigma] = 2.75; M = 25; maxN = 23;

(* Expansion coefficients, 1st set *)

a0[m_] := a0[m] = ((Sqrt[Pi]*(m - 1/2))/(2*M^2*H))*

  Sum[E^(\[Stigma]^2/4 - n^2*H^2)*Sin[(Pi*(m - 1/2)*

    (n*H + \[Stigma]/2))/(M*H)], {n, -maxN, maxN}];

b0[m_] := b0[m] = (-(I/(M*Sqrt[Pi])))*

  Sum[E^(\[Stigma]^2/4 - n^2*H^2)*Cos[(Pi*(m - 1/2)*

    (n*H + \[Stigma]/2))/(M*H)], {n, -maxN, maxN}];

c0[m_] := c0[m] = (Pi*(m - 1/2))/(2*M*H);

(* Equation (7) from Ref. [29] *)

\[CapitalOmega][z_] := \[CapitalOmega][z] = Sum[(a0[m] + b0[m]*z)/

  (c0[m]^2 - z^2), {m, 1, M - 2}];

w1[z_] := \[CapitalOmega][z + I*(\[Stigma]/2)];

The second cell is required to instantiate Equation (8) from [29].

Clear[a1, b1, c1, d1, w2];

(* Expansion coefficients, 2nd set *)

a1[m_] := a1[m] = b0[m]*(((Pi*(m - 1/2))/(2*M*H))^2 -

  (\[Stigma]/2)^2) + I*a0[m]*\[Stigma];

b1[m_] := b1[m] = b0[m];

c1[m_] := c1[m] = (((Pi*(m - 1/2))/(2*M*H))^2 +

  (\[Stigma]/2)^2)^2;

d1[m_] := d1[m] = 2*((Pi*(m - 1/2))/(2*M*H))^2 -

  \[Stigma]^2/2;

(* Equation (8) from Ref. [29] *)

w2[z_] := E^(-z^2) + z*Sum[(a1[m] - b1[m]*z^2)/

  (c1[m] - d1[m]*z^2 + z^4), {m, 1, M - 2}];

This third cell is required to define Equation (9) from [29].

Clear[w3];

(* Equation (9) from Ref. [29] *)

w3[z_] := I/(Sqrt[Pi]*(z - 1/(2*(z - 1/(z - 3/

  (2*(z - 2/(z - 5/(2*(z - 3/(z - 7/(2*(z - 4/

    (z - 9/(2*(z - 5/(z - 11/(2*z))))))))))))))))));

Once these three equations are instantiated, we need to distribute them correspondingly on the complex plane. The forth cell below contains code to accomplish this task.

Clear[wUp, w];

(* Complex error function for upper complex plane *)

wUp[z_] := If[Abs[z] > 8, w3[z],

  If[Im[z] > 0.05*Abs[Re[z]], w1[z], w2[z]]];

(* Complex error function for entire complex plane *)

w[z_] := If[Im[z] >= 0, wUp[z],

  Conjugate[2*E^-Conjugate[z]^2 - wUp[Conjugate[z]]]];

Now the code for computation of the complex error function is ready to use.

It should be noted that another Mathematica code for high-accuracy computation of the complex error function can be downloaded from [33]. This code was written by Jan Mangaldan on the bases of three rational approximations described in our publication [30].

4.2. Fourier Transform

The Mathematica codes below consist of six cells. The code shown in the first cell below defines the Voigt function $K(x,y)$ and imaginary Voigt function $L(x,y)$, the rectangular function $f_r\left(t\right)$ and the sawtooth function $f_s\left(t\right)$.

Clear[K,L,fr,fs];

(* Defining K(x,y) and L(x,y) functions *)

K[x_, y_] := Re[w[x + I*y]];

L[x_, y_] := Im[w[x + I*y]];

(* Rectangular function *)

fr[t_] := 1/((2*t)^(2*35) + 1);

(* Sawtooth function *)

fs[t_] := t*fr[t];

The second cell below generates two look-up tables for the values $V_k(\pi \nu c,nh/c)$ and $V_{\ell }(\pi \nu c,nh/c)$. It also generates a list of grid-points for the parameter $\nu$.

Clear[lookUpTab1, lookUpTab2, nuList];

(* Parameters for computation *)

h = 0.02; c = 0.025; nMax = 25;

(* Computing two look-up tables *)

lookUpTab1 = Table[{K[Pi*\[Nu]*c, (n*h)/c] +

  K[Pi*\[Nu]*c, ((-n)*h)/c]}, {n, 1, nMax},

    {\[Nu], -2*Pi, 2*Pi, 0.1}];

lookUpTab2 = Table[{L[Pi*\[Nu]*c, (n*h)/c] -

  L[Pi*\[Nu]*c, ((-n)*h)/c]}, {n, 1, nMax},

    {\[Nu], -2*Pi, 2*Pi, 0.1}];

nuList = Table[\[Nu], {\[Nu], -2*Pi, 2*Pi, 0.1}];

The code in the following third cell is needed to format two look-up table with values $V_k(\pi \nu c,nh/c)$ and $V_{\ell }(\pi \nu c,nh/c)$ for plotting the graphs. It is also required to join $V_k(\pi \nu c,nh/c)$ and $V_{\ell }(\pi \nu c,nh/c)$ values together with corresponding values of the parameter $\nu$.

Clear[ftList1, ftList2];

(* Main computations by using look up tables *)

ftList1 = Flatten[h*(fr[0]/E^(Pi*nuList*c)^2 +

  Sum[(fr[n*h]*lookUpTab1[[n]])/E^((n*h)/c)^2, {n, 1, nMax}])];

ftList2 = Flatten[h*Sum[(fs[n*h]*lookUpTab2[[n]])/E^((n*h)/c)^2,

  {n, 1, nMax}]];

(* Arranging FT data lists *)

ftList1 = Table[{nuList[[n]], ftList1[[n]]}, {n, 1, Length[nuList]}];

ftList2 = Table[{nuList[[n]], ftList2[[n]]}, {n, 1, Length[nuList]}];

The code in next forth cell is required to generate the references in accordance with equations (35) and (36).

Clear[ftRef1, ftRef2];

(* FT references *)

ftRef1 = Table[{\[Nu], Sinc[Pi*\[Nu]]}, {\[Nu], -2*Pi, 2*Pi, 0.1}];

ftRef2 = Table[{\[Nu], (Pi*\[Nu]*Cos[Pi*\[Nu]] - Sin[Pi*\[Nu]])/(2*

  Pi^2*\[Nu]^2)}, {\[Nu], -2*Pi, 2*Pi, 0.1}];

The code in the fifth cell below produces the graph shown in the Fig. 3.

(* Plotting the graphs from the data lists *)

ListPlot[{ftList1, ftRef1, ftList2, ftRef2}, PlotRange -> All,

  Joined -> True, PlotStyle -> {{Lighter[Green, 0], Thickness[0.005]},

    {Black, Dashed, Thickness[0.0025]}, {Lighter[Magenta, 0.5],

Thickness[0.005]}, {Black, Dashed, Thickness[0.0025]}},

  PlotRange -> {{-2*Pi, 2*Pi}, {-0.3, 1.1}},

    AxesLabel -> {"\[Nu]", None}]

Lastly, the code in following sixth cell applies the derived formula (44) to generate the same Fig. 3 without Voigt functions.

(* Plotting the graphs by using the FT formula (42) *)

f[t_] := fr[t] + fs[t];

ft[\[Nu]_] := (h*(f[0] + Sum[f[n*h]/E^(2*Pi*I*\[Nu]*n*h) +

  f[(-n)*h]*E^(2*Pi*I*\[Nu]*n*h), {n, 1, nMax}]))/

    E^(Pi*\[Nu]*c)^2;

Plot[{Re[ft[\[Nu]]], Sinc[Pi*\[Nu]], Im[ft[\[Nu]]],

  (Pi*\[Nu]*Cos[Pi*\[Nu]] - Sin[Pi*\[Nu]])/(2*(Pi*\[Nu])^2)},

    {\[Nu], -2*Pi, 2*Pi}, PlotRange -> All, PlotStyle ->

      {{Lighter[Green, 0], Thickness[0.005]}, {Black, Dashed,

Thickness[0.0025]}, {Lighter[Magenta, 0.5],

  Thickness[0.005]},{Black, Dashed, Thickness[0.0025]}},

    PlotRange -> {{-2*Pi, 2*Pi}, {-0.3, 1.1}},

      AxesLabel -> {"\[Nu]", None}]

The codes shown in this section can be copy-pasted directly to the Mathematica notebook.

5. Conclusions

An alternative method of rational approximation of the FT based on the real and imaginary parts of the complex error function (12) is developed. Unlike the rational approximation (11) of the FT, the expansion coefficients $V_k(\pi \nu c,nh/c)$ and $V_{\ell }(\pi \nu c,nh/c)$ in this method do not depend on values of the sampled function $f\left(t\right)$. Since the values of the Voigt functions remain always the same, this approach can be used for rapid computation with help of look-up tables. We also show that this rational approximation of the FT can also be rearranged in a trigonometric form (44) with an exponential multiplier $e^{-{\left(\pi \nu c\right)}^2}$ acting like the Hamming window that removes periodicity.