Open access peer-reviewed chapter - ONLINE FIRST

# A Fast Method for Numerical Realization of Fourier Tools

By Anry Nersessian

Submitted: June 16th 2020Reviewed: September 24th 2020Published: December 2nd 2020

DOI: 10.5772/intechopen.94186

## 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.

### Keywords

• Fourier series
• Fourier transforms
• approximation
• over-convergence phenomenon
• numerical methods

## 1. Introduction

One of the classical tools of Mathematics is the apparatus of Fourier series, based on the orthogonal system eiπkx,k=0,±1,±2,complete in L211. 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 11), an intense oscillation arises in their neighborhood, and the uniform convergence is absent (the Gibbs phenomenon). This leads to a slow L2-convergence on the entire segment 11. 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 fis given on the segment 11with Fourier coefficients fs,s=0,±1,,±n,n1,and with the following jumps of the function for its derivatives until degrees of the order q1at the points ak, 1<a1<<am=1,1m<,

Ap,kf=fkap0fkap+0,k=0,1,,q0,p=1,,m.E1

In the neighborhoods of other points we assume that fCq+1. Let us construct a function g=gx,x11with Fourier coefficients gs,0sn, which has the same jumps at the same points, and gCq+1at 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=fgis qtimes continuously differentiable on whole axis, and then

fxgx=s=nnfsgseiπsx+rsx,

where rsx=osq,s,x11.

Therefore, taking into account only first 2n+1Fourier coefficients and truncating the remainder term rs, it is possible to approximate fin the form

fxfnx=defgx+k=nnfsgseiπsxE2

However, the computing of the jumps Askfdirectly by the function fsharply 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 fs(see also [5]). It is easy to obtain using the integration in parts the following asymptotic representation of the Fourier coefficients:

fs=12p=1meiπsapk=0q1Ap,kiπsk+1+rs,rn=osq,s.E3

As a function gfrom (2) K. Eckhoff used those Bernoulli polynomials Bkx,x11,k0, which Fourier coefficients bk,shave the following simple form

bk,s=0,s=0,1s+12iπsk+1,s0,k=1,2,.E4

Denoting B0x=1, polynomials Bkx,k=0,1,,n,x11compose 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

Fnx=defp=1mk=0qAp,kBkxap1+s=nnωseiπsx,E5

where the quantities ωsare given explicitly, converges to fwith the rate onq,n.

K. Eckhoff suggests to find approximate values of jumps A˜p,kAp,kby solving the following system of linear equations with the Vandermonde matrix, obtained by principal part (1) choosing the indexes s=sk,k=1,2,,mq+1,θnskn,0<θ=const<1

fs=12p=1meiπsapk=0qA˜p,kiπsk+1,s=s1,s2,,smq+1E6

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

Snx=defkDnfkexpiπkx,x11,E7

where

fs=1211fxeiπsxdx,s=0,±1,±2,E8

are Fourier coefficients of f, and Dn=dk,k=1,,n,is a set of ndifferent integers (n1). We will assume that D0=Ø.

Definition 2. Letn1be a fixed integer. Consider a system of functionsUn=expiπλkx,λkC,x11,k=1,2,,n, whereλkare arbitrary parameters. Consider the linear spanQn=spanUn. We call a functionqQnas quasi-polynomial of degree at most n.

It is easy to see that qQn, if and only if when either qx0or qx=kPβkxexpiπλkx, where polynomials Pβkx0have the exact degree βkand m=k1+βkn. The number mwill be considered below as the degree of the quasi-polynomial q.

Definition 3. We call the systemλkCas parameters of the quasi-polynomialqx=kPβkxexpiπλkxQnand the number1+βkas multiplicity of the parameterλk.

Definition 4. Let parametersλk,k=1,,nare1fixed. We denote byQnλkthe set of all correspondingqxQn.

Remark 1. The set of quasi-polynomialsqxQnis invariant with respect to a linear change of the variablex. The set of quasi-polynomialsqxQnλkis invariant with respect to a shift of the variablex.

For a fixed integer n1, consider a set of non-integer parameters λkC,kDn, and the following infinite sequence

tr,s=def1srpDnprsprpkDnrλksλk,rDn,s=0,±1,.E9

Remark 2. The last product in(3)is invariant with respect to the numbering of the parametersλk. Here this numbering is tied to the setkDn.

In the general case, the parametersλk(i.e., andtr,s) may depend onn. In order not to complicate the notation, we will indicate this dependence as needed.

Further we denote (see [17] for details)

Trx=defexpiπrx+sDntr,sexpiπsx,rDn,x11,fxFnx=defrDnfrTrx,Rnx=deffxFnxE10

The system Trx1/2expiπrx,rDn, is biorthogonal on the segment x11and L2-error of approximation fxFnxcan be found from the formula

Rn2=sDnfsrDnfrtr,s2E11

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 (x11)

expiπλx=s=sincπsλexpiπsx,λC,E12

where sincz=sinz/z,sinc0=1,zC, plays a key role in the future. This Fourier series can be repeatedly differentiated by the parameter λ.

### 2.2 The case of integer or multiple parameters

To include the case of a quasi-polynomial qQncontaining some integer parameters from λk, consider the following quasi-polynomials associated with formula (12) and sequence (9)

Λj,kx=defi2k2πk1k1!dk1dλjk1cscπλjexpiπλjx=s=1sexpiπsx2sλpk,λjC,jDn,k1,x11.E13

Lemma 1. Ifλjis an integer parameter, then

Λj,kx=BkxexpiπλjxE14

where Bkxis a Bernoulli polynomial (see (4)).

Since the Bernoulli polynomial Bkxcorresponds to one k-multiple parameter λj=0in the KE - method, formula (14) can be considered as a generalization of Bernoulli polynomials to the case of any integer λj.

Remark 3. Like the Bernoulli polynomials in the KE- method, quasi-polynomialsΛj,k,k1can be decomposed into a Fourier series with coefficients(8), for whichfs=Osk,s.

If in the sequence (9) parameters λkand its corresponding multiplicities nkare known, then we will use the representation s=0±1

tr,s=1srpDnprsprpjDmrλjsλjnj,9E15

where rDn,DmDn,nqare corresponding positive integers, and jDmnj=n,λpλqif pq. However, this sequence is still defined only for non-integer numbers in λp.

Consider now the possibility of including in the consideration of some integer parameters in (10). First of all, note that if for some jDn,λjis integer, then tr,λj=. So here the system Trfrom (10) does not exist.

If λjare integers for j,jDn, then, firstly, it is natural to accept tr,r=1, and secondly, we notice that for srthe number of the products in tr,sare reduced.

Finally, we note that in the latter case, the sequence tr,scan be represented as a sum of simple fractions with respect to sC. Bearing in mind Lemma 1, it is not difficult to proof, that.

Lemma 2. Let in (9) the parametersλj,jDm, be given and the subsetΛconsists of all integer parameters fromλj. Then.

1. If there is anλkΛ,λkDn, then the systemTrfrom(10)does not exist.

2. IfΛ=Λ1Λ2Dn, whereΛ1Dmcontains only different integers andΛ2Dmcontains only multiple integers, then

Trx=expiπrx,rΛ1.E16

With aλrΛ2we can apply Lemma 1 and also find the explicit form of functionsTrx.

• IfΛØandΛDnthen the systemTrx,rDm\Λis determined based on the sequence, similar to(15), in whichnj=1,j, andnis replaced bynpwherepis the number of elementsΛ.

• ### 2.3 Explicit form of the system Tr

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 λk.

Theorem 1. Suppose the sequence(15)is given and the possible integer parameters inλksatisfy the conditionλkDm. Then the corresponding functionsTrare quasi-polynomials and have the following explicit form

Trx=jDmk=1njcr,j,kΛj,kx,rDn,x11E17

where

cr,j,k=1r+1iπkpDmrλpnp2njk!pDnprrpdnjkdλjnjkpDnprλjppDnpjλjλpnp

### 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 AlgorithmA.

Let Snxbe a truncated Fourier series (see (7)). Consider formulas (9)(10) with symbolic (not numerical) parameters λq,qDn. Let us now choose a new set of integer numbers D˜n=d˜k,k=1,,n,D˜nDn=Ø. To determine the parameters λq,qDn, we additionally use Fourier coefficients fs,sD˜n, and solve the following system of equations (compare with (11))

fsrDnfrtr,s=0,sD˜nE18

regarding parameters λq. Note that Eq. (18) is essentially nonlinear. If the solution exists, then the system Trfrom (4) is used to approximate fby Fn(see (10)).

Here we show how Algorithm Acan be step-by-step realized using given 2nFourier coefficients fdkfd˜k(see (7)).

Step 1. Using the representation (9), we formally bring the left side of (18) to a common denominator, and fix the conditions λqd˜k,k=1,2,,n. Then Eq. (18) will take the form (sD˜n)

fsqDnsλq=rDn1sdrfrsincdrλrpDnpdrsdpdrdpqDndrλq.E19

Step 2. Using the Vieta’s formula decompose the products in (18) containing the parameters λkand denote

ek=1ki1<i2<<ikλi1λi2,,λik,k=1,2,,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 ek.

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

Step 4. Again, according to the Vieta’s formulae find all rootszk,k=1,2,,nof the corresponding polynomial Pz=p=0nepznp,e0=1, and put λk=zk.

Step 5. Realize the approximation fFn(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}. LetfQn(see Definition 2) and the setsDn,D˜nand the Fourier coefficientsfs,sDnD˜n, of thefbe given. Denote byΛthe set of integer parameters in the representationf=kPmkxexpiπμkx. In order for the approximation by AlgorithmAto be exact (that isfxFnx,x11), it is necessary and sufficient thatΛDnD˜n.

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

Definition 5. LetfL211and the setsDn,D˜nand the Fourier coefficientsfs,sDnD˜n, of thefbe given. Applying Algorithm 1, we get parametersλk,kDn. The accelerating the convergence of truncated Fourier series

S2nx=kDnD˜nfkexpiπkx,x11,E20

we define by formula (17).

Remark 5. There are no Fourier coefficientsfr,rD˜minformula (17). These coefficients are used in AlgorithmAto an optimal choice of parametersλk,kDm.

## 3. An numerical algorithm for Fourier transform

Consider the Fourier transform Fof a function fL2the form

Fω=12πfteitωdt,ωE21

### 3.1 Algorithm construction scheme

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

Fω=12πk=khk+hfteitωdtE22

On the other hand, our method allows one to approximate function ffrTron each segment khk+h,k,using Algorithm Aand 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 Trxis a linear combination of functions of the form expiπλ1xor its derivative with respect to λof a certain order. For example, if k=0, then for functions expiπλ1x, we have (see (12), here k=0)

hhexpiπλ1teitωdt=2hsinchπλ1+ω

On the complex plane sC, 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 functionfinformula (21), intervalcan 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 n1,Dn=01n1,D˜n=nn+11.

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 (see Remark 6).

3. The computer errs or freezes most often due to type 0/0uncertainty. 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 Awe have λ1=1.00012,λ2=1.00031. At this point, so that the calculations do not stop, it is necessary to make following correction in the code: λ1=λ2=1. Of course, here we mean the permissible error 103when computer precision is 105(see above item 1).

1. During the operation of the proposed algorithm for Fourier transform, L2- error can be monitored every time the adaptive algorithm is applied. In practice, it has the same order as the L2-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)).

fx=e4x13+i2e4x122,xE23
Fω=42πe3ω2+16+ie116ωω8i42.ωE24

Our goal is to provide a final value of a L2- error that does not exceed 102. With the help of numerical integration (in this case, explicitly), one can make sure that

32fx2dx+2fx2dx=4.×103

and when deleting these semi-infinite intervals, further actions should lead to an error no greater than 5.99×103.

On the other hand, we see that function fis 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: I1=21/3,I2=1/33/2. Method of item 3.1 at n=3was applied to each of these intervals.

Monitoring showed (see item 4) that L2- error for I1is 1.×103, and for I2is 1.44×103. We see that this is less than allowed above.

The total L2- error turned out to be 4.38×103. As for the uniform error, it is reached approximately at point ω=0.185and is equal to 1.04×103.

These calculations were performed using system Wolfram Mathematica 9.

## 4. Conclusion

Those who have mastered the technique of using Algorithm Acan easily obtain a similar method of sine and cosine Fourier transforms. Despite the fact that with an increase in nthe 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]).

## Notes

• Unless otherwise stated, we will not exclude the repetition of values of λk in a set λk.

chapter PDF

## More

© 2020 The Author(s). Licensee IntechOpen. This chapter is distributed under the terms of the Creative Commons Attribution 3.0 License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

## How to cite and reference

### Cite this chapter Copy to clipboard

Anry Nersessian (December 2nd 2020). A Fast Method for Numerical Realization of Fourier Tools [Online First], IntechOpen, DOI: 10.5772/intechopen.94186. Available from: