A Fast Method for Numerical Realization of Fourier Tools

Abstract

This chapter presents new application of author’s recent algorithms for fast summations of truncated Fourier series. A complete description of this method is given, and an algorithm for numerical implementation with a given accuracy for the Fourier transform is proposed.

1. Introduction

One of the classical tools of Mathematics is the apparatus of Fourier series, based on the orthogonal system $\left\{e^{\mathit{\text{iπkx}}}\right\},k=0,\pm 1,\pm 2,\dots$ complete in $L_2\left[-1,1\right]$. However, in practice it is extremely limited because of the poor approximation of piecewise smooth functions. Thus, in the case when a function has discontinuity points (taking into account also discontinuities at the ends of the interval $\left[-1,1\right]$), an intense oscillation arises in their neighborhood, and the uniform convergence is absent (the Gibbs phenomenon). This leads to a slow $L_2$-convergence on the entire segment $\left[-1,1\right]$. The corresponding phenomenon is observed in the case of Fourier interpolation. Classical methods of summation truncated Fourier series (see, for example, [1], Chapter 3) do not actually change the situation.

Since the beginning of the last century, a huge amount of research has been carried out, devoted to the methods of effective summation of the Fourier series. Let us briefly dwell on the works to overcome this situation (for detiles see, for example, articles [2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13] and their links).

The pioneer of “overcoming the Gibbs phenomenon” is A.N.Krylov, who at the dawn of the 20th century (see [2]) proposed methods that he later developed in the monograph [3]. In particular, he proposed the following approach. Let a piecewise smooth function $f$ is given on the segment $\left[-1,1\right]$ with Fourier coefficients $\left\{f_s\right\},s=0,\pm 1,\dots ,\pm n,n\ge 1,$ and with the following jumps of the function $f$ or its derivatives until degrees of the order $q\ge 1$ at the points $\left\{a_k\right\}$, $-1<a_1<\cdots <a_m=1,1\le m<\infty ,$

$$\begin{array}{c}{A}_{p,k}\left(f\right)=f^{\left(k\right)}\left(a_p-0\right)-f^{\left(k\right)}\left(a_p+0\right),\\ k=0,1,\cdots ,q\ge 0,p=1,\cdots ,m.\end{array}$$E1

In the neighborhoods of other points we assume that $f\in C^{q+1}$. Let us construct a function $g=g\left(x\right),x\in \left[-1,1\right]$ with Fourier coefficients $\left\{g_s\right\},0\le \mid s\mid \le n$, which has the same jumps at the same points, and $g\in C^{q+1}$ at the neighborhoods of other points. Known the jumps (1), one can construct, e.g. piecewise-polynomial $g$. As a result, 2-periodic extension of the function $F=\left(f-g\right)$ is $q$ times continuously differentiable on whole axis, and then

$$f\left(x\right)-g\left(x\right)=\sum _{s=-n}^n\left(f_s-g_s\right)e^{\mathit{\text{iπsx}}}+r_s\left(x\right),$$

where $r_s\left(x\right)=o\left(s^{-q}\right),s\to \infty ,x\in \left[-1,1\right]$.

Therefore, taking into account only first $\left(2n+1\right)$ Fourier coefficients and truncating the remainder term $r_s$, it is possible to approximate $f$ in the form

$$f\left(x\right)\simeq f_n\left(x\right)\stackrel{\mathit{def}}{=}g\left(x\right)+\sum _{k=-n}^n\left(f_s-g_s\right)e^{\mathit{\text{iπsx}}}$$E2

However, the computing of the jumps $\left\{A_{\mathit{sk}}\left(f\right)\right\}$ directly by the function $f$ sharply limits the scope of practical application of the method of A. N. Krylov.

The “spectral” method proposed by Kurt Eckhoff in [4] (1993) turned out to be more practical, as it is based only on the use of coefficients $\left\{f_s\right\}$ (see also [5]). It is easy to obtain using the integration in parts the following asymptotic representation of the Fourier coefficients:

$$f_s=-\frac{1}{2}\sum _{p=1}^m{e}^{-\mathit{i\pi s}{a}_p}\sum _{k=0}^{q-1}\frac{A_{p,k}}{{\left(\mathit{i\pi s}\right)}^{k+1}}+r_s,r_n=o\left(s^{-q}\right),s\to \infty .$$E3

As a function $g$ from (2) K. Eckhoff used those Bernoulli polynomials $\left\{B_k\left(x\right)\right\},x\in \left[-1,1\right],k\ge 0$, which Fourier coefficients $\left\{b_{k,s}\right\}$ have the following simple form

$$b_{k,s}=\left\{\begin{array}{ll}0,& s=0,\\ \frac{{\left(-1\right)}^{s+1}}{2{\left(\mathit{i\pi s}\right)}^{k+1}},& s\ne 0,k=1,2,.\dots \end{array}\right.$$E4

Denoting $B_0\left(x\right)=1$, polynomials $\left\{B_k\left(x\right)\right\},k=0,1,\dots ,n,x\in \left[-1,1\right]$ compose a basis on the space of polynomials of degree $n$. Bernoulli polynomials extended to the real axis with period 2 as piecewise-smooth functions.

According to Krylov’s scheme, the sequence

$$F_n\left(x\right)\stackrel{\mathit{def}}{=}\sum _{p=1}^m\sum _{k=0}^q{A}_{p,k}{B}_k\left(x-a_p-1\right)+\sum _{s=-n}^n{\omega }_s{e}^{\mathit{\text{iπsx}}},$$E5

where the quantities $\left\{{\omega }_s\right\}$ are given explicitly, converges to $f$ with the rate $o\left(n^{-q}\right),n\to \infty$.

K. Eckhoff suggests to find approximate values of jumps $\left\{{\tilde{A}}_{p,k}\approx A_{p,k}\right\}$ by solving the following system of linear equations with the Vandermonde matrix, obtained by principal part (1) choosing the indexes $s=s_k,k=1,2,\dots ,m\left(q+1\right),\mathit{\theta n}\le \mid s_k\mid \le n,0<\theta =\mathit{\text{const}}<1$

$$f_s=-\frac{1}{2}\sum _{p=1}^m{e}^{-i\mathit{\pi s}{a}_p}\sum _{k=0}^q\frac{{\tilde{A}}_{p,k}}{{\left(\mathit{i\pi s}\right)}^{k+1}},s=s_1,s_2,\dots ,s_{m\left(q+1\right)}$$E6

Further, the problem is solved by the Krylov method. It is natural to call this acceleration scheme the Krylov-Eckhoff method (KE-method). Over the past two to three decades, in the KE-method scheme, not only polynomials have been used as functions of $g$.

In [14], for a smooth function $f$, a combination of the solution of Eq.(6) and the Pade approximation (applied to the asymptotic expansion of the Fourier coefficients) is implemented there. As a result of numerical experiments it turned out that, as a rule, such an algorithm is “almost exact” on certain infinite-dimensional spaces, although only a finite number of Fourier coefficients is used. In [15], a similar phenomenon was discovered for Bessel series, and in [16] for some eigenfunction expansions. The theoretical substantiation of this phenomenon for Fourier series and the corresponding algorithm were published in [17, 18, 19]. Our goal is to use this method in the case of one-dimensional Fourier transforms. Let us start with a complete description of the required algorithms from [17].

2. Acceleration of convergence of Fourier series

2.1 Basic definitions

The classical definition of the partial sums of the Fourier series is based on a gradual increase in the frequencies of the Fourier system. Below we use a more general notations from the work [17].

Definition 1. Let us call the truncated Fourier series any sum of the form

$$S_n\left(x\right)\stackrel{\mathit{def}}{=}\sum _{k\in D_n}{f}_{k}exp\left(i\pi kx\right),x\in \left[-1,1\right],$$E7

where

$$f_s=\frac{1}{2}{\int }_{-1}^{1}f\left(x\right)e^{-i\pi sx}\mathit{dx},s=0,\pm 1,\pm 2,\dots$$E8

are Fourier coefficients of $f$, and $D_n=\left\{d_k\right\},k=1,\dots ,n,$ is a set of $n$ different integers ($n\ge 1$). We will assume that $D_0=\mathrm{Ø}$.

Definition 2. Let$n\ge 1$be a fixed integer. Consider a system of functions$U_n=\left\{exp\left(i\pi {\lambda }_{k}x\right)\right\},{\lambda }_k\in \mathbb{C},x\in \left[-1,1\right],k=1,2,\dots ,n$, where$\left\{{\lambda }_k\right\}$are arbitrary parameters. Consider the linear span$Q_n=\mathit{\text{span}}\left\{U_n\right\}$. We call a function$q\in Q_n$as quasi-polynomial of degree at most n.

It is easy to see that $q\in Q_n$, if and only if when either $q\left(x\right)\equiv 0$ or $q\left(x\right)={\sum }_k{P}_{{\beta }_k}\left(x\right)exp\left(i\pi {\lambda }_{k}x\right)$, where polynomials $P_{{\beta }_k}\left(x\right)\not\equiv 0$ have the exact degree ${\beta }_k$ and $m={\sum }_k\left(1+{\beta }_k\right)\le n$. The number $m$ will be considered below as the degree of the quasi-polynomial $q$.

Definition 3. We call the system$\left\{{\lambda }_k\right\}\subset \mathbb{C}$as parameters of the quasi-polynomial$q\left(x\right)={\sum }_k{P}_{{\beta }_k}\left(x\right)exp\left(i\pi {\lambda }_{k}x\right)\in Q_n$and the number$1+{\beta }_k$as multiplicity of the parameter${\lambda }_k$.

Definition 4. Let parameters$\left\{{\lambda }_k\right\},k=1,\dots ,n$are¹fixed. We denote by$Q_n\left(\left\{{\lambda }_k\right\}\right)$the set of all corresponding$q\left(x\right)\in Q_n$.

Remark 1. The set of quasi-polynomials$q\left(x\right)\in Q_n$is invariant with respect to a linear change of the variable$x$. The set of quasi-polynomials$q\left(x\right)\in Q_n\left(\left\{{\lambda }_k\right\}\right)$is invariant with respect to a shift of the variable$x$.

For a fixed integer $n\ge 1$, consider a set of non-integer parameters $\left\{{\lambda }_k\right\}\subset \mathbb{C},k\in D_n$, and the following infinite sequence

$$t_{r,s}\stackrel{\mathit{def}}{=}{\left(-1\right)}^{s-r}\left(\prod _{\begin{array}{c}p\in D_n\\ p\ne r\end{array}}\frac{s-p}{r-p}\right)\prod _{k\in D_n}\frac{r-{\lambda }_k}{s-{\lambda }_k},r\in D_n,s=0,\pm 1,\dots .$$E9

Remark 2. The last product in(3)is invariant with respect to the numbering of the parameters$\left\{{\lambda }_k\right\}$. Here this numbering is tied to the set$k\in D_n$.

In the general case, the parameters$\left\{{\lambda }_k\right\}$(i.e., and$t_{r,s}$) may depend on$n$. In order not to complicate the notation, we will indicate this dependence as needed.

Further we denote (see [17] for details)

$$\begin{array}{l}{T}_r\left(x\right)\stackrel{\mathit{def}}{=}exp\left(i\pi rx\right)+\sum _{s\cancel{\in }{D}_n}{t}_{r,s}exp\left(i\pi sx\right),r\in D_n,x\in \left[-1,1\right],\\ f\left(x\right)\simeq F_n\left(x\right)\stackrel{\mathit{def}}{=}\sum _{r\in D_n}{f}_r{T}_r\left(x\right),R_n\left(x\right)\stackrel{\mathit{def}}{=}f\left(x\right)-F_n\left(x\right)\end{array}$$E10

The system $\left\{T_r\left(x\right),1/2exp\left(i\pi rx\right)\right\},r\in D_n$, is biorthogonal on the segment $x\in \left[-1,1\right]$ and $L_2$-error of approximation $f\left(x\right)\simeq F_n\left(x\right)$ can be found from the formula

$${‖R_n‖}^2=\sum _{s\cancel{\in }{D}_n}{‖f_s-\sum _{r\in D_n}{f}_r{t}_{r,s}‖}^2$$E11

Here we confine ourselves to the one-dimensional case. A similar universal algorithm for multivariate truncated Fourier series was proposed and numerically implemented in [1].

The following obvious formula ($x\in \left[-1,1\right]$)

$$exp\left(i\pi \lambda x\right)=\sum _{s=-\infty }^{\infty }sinc\left(\pi \left(s-\lambda \right)\right)exp\left(i\pi sx\right),\lambda \in \mathbb{C},$$E12

where $sinc\left(z\right)=sin\left(z\right)/z,sinc\left(0\right)=1,z\in \mathbb{C}$, plays a key role in the future. This Fourier series can be repeatedly differentiated by the parameter $\lambda$.

2.2 The case of integer or multiple parameters

To include the case of a quasi-polynomial $q\in Q_n$ containing some integer parameters from $\left\{{\lambda }_k\right\}$, consider the following quasi-polynomials associated with formula (12) and sequence (9)

$$\begin{array}{l}{\mathrm{\Lambda }}_{j,k}\left(x\right)\stackrel{\mathit{def}}{=}\\ \frac{i^{2-k}}{2{\pi }^{k-1}\left(k-1\right)!}\frac{d^{k-1}}{d{\lambda }_j^{k-1}}\left(csc\left({\mathit{\pi \lambda }}_j\right)exp\left(i\pi {\lambda }_{j}x\right)\right)=\\ \sum _{s=-\infty }^{\infty }\frac{{\left(-1\right)}^{s}exp\left(i\pi sx\right)}{2{\left(\mathit{i\pi }\left(s-{\lambda }_p\right)\right)}^k},\\ {\lambda }_j\in \mathbb{C},j\in D_n,k\ge 1,x\in \left[-1,1\right].\end{array}$$E13

Lemma 1. If${\lambda }_j$is an integer parameter, then

$${\mathrm{\Lambda }}_{j,k}\left(x\right)=B_k\left(x\right)exp\left(i\pi {\lambda }_{j}x\right)$$E14

where $B_k\left(x\right)$ is a Bernoulli polynomial (see (4)).

Since the Bernoulli polynomial $B_k\left(x\right)$ corresponds to one $k$ -multiple parameter ${\lambda }_j=0$ in the KE - method, formula (14) can be considered as a generalization of Bernoulli polynomials to the case of any integer ${\lambda }_j$ .

Remark 3. Like the Bernoulli polynomials in the KE- method, quasi-polynomials$\left\{{\mathrm{\Lambda }}_{j,k}\right\},k\ge 1$can be decomposed into a Fourier series with coefficients(8), for which$f_s=O\left({\left|s\right|}^{-k}\right),s\to \infty$.

If in the sequence (9) parameters $\left\{{\lambda }_k\right\}$ and its corresponding multiplicities $\left\{n_k\right\}$ are known, then we will use the representation $\left(s=0,\pm 1,\dots \right)$

$$t_{r,s}={\left(-1\right)}^{s-r}\left(\prod _{\begin{array}{c}p\in D_n\\ p\ne r\end{array}}\frac{s-p}{r-p}\right)\prod _{j\in D_m}{\left(\frac{r-{\lambda }_j}{s-{\lambda }_j}\right)}^{n_j},\left(9^{\ast }\right)$$E15

where $r\in D_n,D_m\subset D_n,\left\{n_q\right\}$ are corresponding positive integers, and ${\sum }_{j\in D_m}{n}_j=n,{\lambda }_p\ne {\lambda }_q$ if $p\ne q$. However, this sequence is still defined only for non-integer numbers in $\left\{{\lambda }_p\right\}$.

Consider now the possibility of including in the consideration of some integer parameters in (10). First of all, note that if for some $j\notin D_n,{\lambda }_j$ is integer, then $t_{r,{\lambda }_j}=\infty$. So here the system $\left\{T_r\right\}$ from (10) does not exist.

If ${\lambda }_j$ are integers for $j,j\in D_n$, then, firstly, it is natural to accept $t_{r,r}=1$, and secondly, we notice that for $s\ne r$ the number of the products in $t_{r,s}$ are reduced.

Finally, we note that in the latter case, the sequence $t_{r,s}$ can be represented as a sum of simple fractions with respect to $s\in \mathbb{C}$. Bearing in mind Lemma 1, it is not difficult to proof, that.

Lemma 2. Let in ($9^{\ast }$) the parameters$\left\{{\lambda }_j\right\},j\in D_m$, be given and the subset$\mathrm{\Lambda }$consists of all integer parameters from$\left\{{\lambda }_j\right\}$. Then.

1. If there is an${\lambda }_k\in \mathrm{\Lambda },{\lambda }_k\notin D_n$, then the system$\left\{T_r\right\}$from(10)does not exist.

2. If$\mathrm{\Lambda }=\left({\mathrm{\Lambda }}_1\cup {\mathrm{\Lambda }}_2\right)\subset D_n$, where${\mathrm{\Lambda }}_1\subset D_m$contains only different integers and${\mathrm{\Lambda }}_2\subset D_m$contains only multiple integers, then

   $$T_r\left(x\right)=exp\left(i\pi rx\right),r\in {\mathrm{\Lambda }}_1.$$E16

   With a${\lambda }_r\in {\mathrm{\Lambda }}_2$we can apply Lemma 1 and also find the explicit form of functions$\left\{T_r\left(x\right)\right\}$.

3. If$\mathrm{\Lambda }\ne \mathrm{Ø}$and$\mathrm{\Lambda }\ne D_n$then the system$\left\{T_r\left(x\right)\right\},r\in D_m\backslash \mathrm{\Lambda }$is determined based on the sequence, similar to(15), in which$n_j=1,\forall j$, and$n$is replaced by$\left(n-p\right)$where$p$is the number of elements$\mathrm{\Lambda }$.

2.3 Explicit form of the system $\left\{T_r\right\}$

The following result is a generalization of Theorem 1 in [17] to the case that in formula (15) there are integer values among the parameters $\left\{{\lambda }_k\right\}$.

Theorem 1. Suppose the sequence(15)is given and the possible integer parameters in$\left\{{\lambda }_k\right\}$satisfy the condition${\lambda }_k\in D_m$. Then the corresponding functions$\left\{T_r\right\}$are quasi-polynomials and have the following explicit form

$$T_r\left(x\right)=\sum _{j\in D_m}\sum _{k=1}^{n_j}{c}_{r,j,k}{\mathrm{\Lambda }}_{j,k}\left(x\right),r\in D_n,x\in \left[-1,1\right]$$E17

where

$$\begin{array}{l}{c}_{r,j,k}=\frac{{\left(-1\right)}^{r+1}{\left(i\pi \right)}^k{\prod }_{p\in D_m}{\left(r-{\lambda }_p\right)}^{n_p}}{2\left(n_j-k\right)!{\prod }_{\begin{array}{c}p\in D_n\\ p\ne r\end{array}}\left(r-p\right)}\\ \frac{d^{n_j-k}}{d{\lambda }_j^{n_j-k}}\left(\frac{{\prod }_{\begin{array}{c}p\in D_n\\ p\ne r\end{array}}\left({\lambda }_j-p\right)}{{\prod }_{\begin{array}{c}p\in D_n\\ p\ne j\end{array}}{\left({\lambda }_j-{\lambda }_p\right)}^{n_p}}\right)\end{array}$$

2.4 The adaptive algorithm and over-convergence phenomenon

In the main Theorem of paper [17], the phenomenon of over-convergence was theoretically substantiated. The basis of this result was the following Algorithm$\mathfrak{A}$.

Let $S_n\left(x\right)$ be a truncated Fourier series (see (7)). Consider formulas (9)–(10) with symbolic (not numerical) parameters $\left\{{\lambda }_q\right\},q\in D_n$. Let us now choose a new set of integer numbers ${\tilde{D}}_n=\left\{{\tilde{d}}_k\right\},k=1,\dots ,n,{\tilde{D}}_n\cap D_n=\mathrm{Ø}$. To determine the parameters $\left\{{\lambda }_q\right\},q\in D_n$, we additionally use Fourier coefficients $\left\{f_s\right\},s\in {\tilde{D}}_n$, and solve the following system of equations (compare with (11))

$$f_s-\sum _{r\in D_n}{f}_r{t}_{r,s}=0,s\in {\tilde{D}}_n$$E18

regarding parameters $\left\{{\lambda }_q\right\}$. Note that Eq. (18) is essentially nonlinear. If the solution exists, then the system $\left\{T_r\right\}$ from (4) is used to approximate $f$ by $F_n$ (see (10)).

Here we show how Algorithm $\mathfrak{A}$ can be step-by-step realized using given $2n$ Fourier coefficients $\left\{f_{d_k}\right\}\bigcup \left\{\left\{f_{{\tilde{d}}_k}\right\}\right\}$ (see (7)).

Step 1. Using the representation (9), we formally bring the left side of (18) to a common denominator, and fix the conditions ${\lambda }_q\ne {\tilde{d}}_k,k=1,2,\dots ,n$. Then Eq. (18) will take the form ($s\in {\tilde{D}}_n$)

$$\begin{array}{l}{f}_s\prod _{q\in D_n}\left(s-{\lambda }_q\right)=\\ \sum _{r\in D_n}\frac{{\left(-1\right)}^{s-d_r}{f}_r}{sinc\left(d_r-{\lambda }_r\right)}\left({\prod }_{\begin{array}{c}p\in D_n\\ p\ne d_r\end{array}}\frac{s-d_p}{d_r-d_p}\right)\prod _{q\in D_n}\left(d_r-{\lambda }_q\right).\end{array}$$E19

Step 2. Using the Vieta’s formula decompose the products in (18) containing the parameters $\left\{{\lambda }_k\right\}$ and denote

$$e_k={\left(-1\right)}^k\sum _{i_1<i_2<\dots <i_k}{\lambda }_{i_1}{\lambda }_{i_2},\dots ,{\lambda }_{i_k},k=1,2,\dots ,n.$$

It is not difficult to see that, as a result, Eq.(18) will take the form of a linear system of equations with respect to the variables $\left\{e_k\right\}$.

Step 3. Solve resulting equation by the least square method.

Step 4. Again, according to the Vieta’s formulae find all roots$\left\{z_k\right\},k=1,2,\dots ,n$ of the corresponding polynomial $P\left(z\right)={\sum }_{p=0}^n{e}_p{z}^{n-p},e_0=1$, and put $\left\{{\lambda }_k\right\}=\left\{z_k\right\}$.

Step 5. Realize the approximation $f\simeq F_n$ (see (10) and (17)).

Remark 4. Least square method provides only one solution in Step 3.

The following theorem is the main result of the paper [17].

Theorem 2. {The phenomenon of the over-convergence}. Let$f\in Q_n$(see Definition 2) and the sets$D_n$,${\tilde{D}}_n$and the Fourier coefficients$\left\{f_s\right\},s\in D_n\cup {\tilde{D}}_n$, of the$f$be given. Denote by$\mathrm{\Lambda }$the set of integer parameters in the representation$f={\sum }_k{P}_{m_k}\left(x\right)exp\left(i\pi {\mu }_{k}x\right)$. In order for the approximation by Algorithm$\mathfrak{A}$to be exact (that is$f\left(x\right)\equiv F_n\left(x\right),x\in \left[-1,1\right]$), it is necessary and sufficient that$\mathrm{\Lambda }\subseteq D_n\cup {\tilde{D}}_n$.

It is now natural to formulate the following definition of the acceleration of the convergence of a truncated Fourier series.

Definition 5. Let$f\in L_2\left[-1,1\right]$and the sets$D_n$,${\tilde{D}}_n$and the Fourier coefficients$\left\{f_s\right\},s\in D_n\cup {\tilde{D}}_n$, of the$f$be given. Applying Algorithm 1, we get parameters$\left\{{\lambda }_k\right\},k\in D_n$. The accelerating the convergence of truncated Fourier series

$$S_{2n}\left(x\right)=\sum _{k\in D_n\cup {\tilde{D}}_n}{f}_{k}exp\left(i\pi kx\right),x\in \left[-1,1\right],$$E20

we define by formula (17).

Remark 5. There are no Fourier coefficients$\left\{f_r\right\},r\in {\tilde{D}}_m$informula (17). These coefficients are used in Algorithm$\mathfrak{A}$to an optimal choice of parameters$\left\{{\lambda }_k\right\},k\in D_m$.

3. An numerical algorithm for Fourier transform

Consider the Fourier transform $F$ of a function $f\in L_2\left(-\infty ,\infty \right)$ the form

$$F\left(\omega \right)=\frac{1}{\sqrt{2\pi }}{\int }_{-\infty }^{\infty }f\left(t\right)e^{it\omega }\mathit{dt},\omega \in \left(-\infty ,\infty \right)$$E21

3.1 Algorithm construction scheme

The method of the item 2 can be applied by linear change of variable to function $f$ defined on any finite segment. For a fixed $h>0$formula (21) can be rewritten as

$$F\left(\omega \right)=\frac{1}{\sqrt{2\pi }}\sum _{k=-\infty }^{\infty }{\int }_{k-h}^{k+h}f\left(t\right)e^{it\omega }\mathit{dt}$$E22

On the other hand, our method allows one to approximate function $f\simeq \sum f_r{T}_r$ on each segment $\left(k-h,k+h\right),\forall k,$ using Algorithm $\mathfrak{A}$ and formula (17).

Our idea is as follows. First, based on the required 1 error, you can leave a fnite number of intervals. Second, on the remaining intervals, the integrals are calculated explicitly. Really, each function $T_r\left(x\right)$ is a linear combination of functions of the form $exp\left(i\pi {\lambda }_{1}x\right)$ or its derivative with respect to $\lambda$ of a certain order. For example, if $k=0$, then for functions $exp\left(i\pi {\lambda }_{1}x\right)$, we have (see (12), here $k=0$)

$${\int }_{-h}^{h}exp\left(i\pi {\lambda }_{1}t\right)e^{it\omega }\mathit{dt}=2hsinc\left(h\left({\mathit{\pi \lambda }}_1+\omega \right)\right)$$

On the complex plane $s\in \mathbb{C}$, this is an entire function of exponential type. It is easy to see that this is true in the general case.

Remark 6. Depending on function$f$informula (21), interval$\left(-\infty ,\infty \right)$can be split into a finite number of unequal intervals instead of partition(22).

Let us turn to the main details of the algorithm for such a numerical implementation of the Fourier transform for a piecewise smooth function $f$.

3.2 Implementation notes

To make our approach understandable to a wide range of specialists, when developing an applied algorithm, we will assume that (see formulas (7–19)) for a fixed $n\ge 1,D_n=\left\{0,1,\dots ,n-1\right\},{\tilde{D}}_n=\left\{-n,-n+1,\dots ,-1\right\}$.

Consider now the following when writing code to implement the one-dimensional Fourier transform.

1. The computer used must be installed with programs that can perform symbolic mathematical operations as well as numerical integration (for example, as system Wolfram Mathematica, see [19]). It is desirable to use the highest possible accuracy of calculations.

2. The singularities of piecewise-smooth function $f$ (discontinuity points of the function and the discontinuity of the derivatives of low order) must be included in the partition of interval $\left(-\infty ,\infty \right)$ (see Remark 6).

3. The computer errs or freezes most often due to type $0/0$ uncertainty. The point is that in practice, due to rounding errors, when the parameters coincide, they are often not fixed (see formula (15)). This must be taken into account also in the presence of integer parameters (see item 2.2).

For example, let in the process of computer operation at the output of Algorithm $\mathfrak{A}$ we have ${\lambda }_1=1.00012,{\lambda }_2=1.00031$. At this point, so that the calculations do not stop, it is necessary to make following correction in the code: ${\lambda }_1={\lambda }_2=1$. Of course, here we mean the permissible error ${10}^{-3}$ when computer precision is ${10}^{-5}$ (see above item 1).

1. During the operation of the proposed algorithm for Fourier transform, $L_2$- error can be monitored every time the adaptive algorithm is applied. In practice, it has the same order as the $L_2$-error for the Fourier transform.

3.3 A simple example

Let us explain the details on the following example of an approximate finding of the Fourier transform (see notation (21)).

$$f\left(x\right)=e^{-4\left|x-\frac{1}{3}\right|}+\frac{i}{2}{e}^{-4{\left(x-\frac{1}{2}\right)}^2},x\in \left(-\infty ,\infty \right)$$E23

$$F\left(\omega \right)=\frac{4\sqrt{\frac{2}{\pi }}{e}^{\frac{\mathit{i\omega }}{3}}}{{\omega }^2+16}+\frac{{\mathit{ie}}^{-\frac{1}{16}\omega \left(\omega -8i\right)}}{4\sqrt{2}}.\omega \in \left(-\infty ,\infty \right)$$E24

Our goal is to provide a final value of a $L_2$ - error that does not exceed ${10}^{-2}$. With the help of numerical integration (in this case, explicitly), one can make sure that

$$\sqrt{{\int }_{\frac{3}{2}}^{\infty }{\left|f\left(x\right)\right|}^2\mathit{dx}+{\int }_{-\infty }^{-2}{\left|f\left(x\right)\right|}^2\mathit{dx}}=4.\times {10}^{-3}$$

and when deleting these semi-infinite intervals, further actions should lead to an error no greater than $5.99\times {10}^{-3}$.

On the other hand, we see that function $f$ is infinitely differentiable everywhere except the point $x=1/3$, where its derivative has a jump. Therefore (see item 2), we can consider the following partition into two intervals: $I_1=\left(-2,1/3\right),I_2=\left(1/3,3/2\right)$. Method of item 3.1 at $n=3$ was applied to each of these intervals.

Monitoring showed (see item 4) that $L_2$- error for $I_1$ is $1.\times {10}^{-3}$, and for $I_2$ is $1.44\times {10}^{-3}$. We see that this is less than allowed above.

The total $L_2$- error turned out to be $4.38\times {10}^{-3}$. As for the uniform error, it is reached approximately at point $\omega =0.185$ and is equal to $1.04\times {10}^{-3}$.

These calculations were performed using system Wolfram Mathematica 9.

4. Conclusion

Those who have mastered the technique of using Algorithm $\mathfrak{A}$ can easily obtain a similar method of sine and cosine Fourier transforms. Despite the fact that with an increase in $n$ the complexity of the corresponding algorithms increases markedly, modern computers allow you to perform the necessary actions with increased bit depth. The proposed method for finding the one-dimensional Fourier transform with a given accuracy can be useful in such different fields of activity as mathematical physics, astronomy, engineering, economics, biology, etc.

There are good prerequisites for other classical (in a broad sense) Fourier tools for efficient applications based on the method of item 2 (see articles [14, 15]).