Open access peer-reviewed chapter

A Fast Method for Numerical Realization of Fourier Tools

Written By

Anry Nersessian

Submitted: 16 June 2020 Reviewed: 24 September 2020 Published: 02 December 2020

DOI: 10.5772/intechopen.94186

From the Edited Volume

Real Perspective of Fourier Transforms and Current Developments in Superconductivity

Edited by Juan Manuel Velazquez Arcos

Chapter metrics overview

663 Chapter Downloads

View Full Metrics


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.


  • 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 f is given on the segment 11 with Fourier coefficients fs,s=0,±1,,±n,n1, and with the following jumps of the function f or its derivatives until degrees of the order q1 at the points ak, 1<a1<<am=1,1m<,


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


where rsx=osq,s,x11.

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


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


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


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


where the quantities ωs are given explicitly, converges to f with the rate onq,n.

K. Eckhoff suggests to find approximate values of jumps A˜p,kAp,k by 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


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




are Fourier coefficients of f, and Dn=dk,k=1,,n, is a set of n different 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 qx0 or qx=kPβkxexpiπλkx, where polynomials Pβkx0 have the exact degree βk and m=k1+βkn. The number m will 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


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)


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


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)


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 qQn containing some integer parameters from λk, consider the following quasi-polynomials associated with formula (12) and sequence (9)


Lemma 1. Ifλjis an integer parameter, then


where Bkx is a Bernoulli polynomial (see (4)).

Since the Bernoulli polynomial Bkx corresponds to one k -multiple parameter λj=0 in 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 λk and its corresponding multiplicities nk are known, then we will use the representation s=0±1


where rDn,DmDn,nq are corresponding positive integers, and jDmnj=n,λpλq if 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,λj is integer, then tr,λj=. So here the system Tr from (10) does not exist.

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

Finally, we note that in the latter case, the sequence tr,s can 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


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

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




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 Snx be 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))


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

Here we show how Algorithm A can be step-by-step realized using given 2n Fourier 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)


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


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,,n of 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


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 F of a function fL2 the form


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>0formula (21) can be rewritten as


On the other hand, our method allows one to approximate function ffrTr on each segment khk+h,k, using Algorithm 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 Trx is a linear combination of functions of the form expiπλ1x or 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)


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/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 A we 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 103 when 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)).


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


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 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: I1=21/3,I2=1/33/2. Method of item 3.1 at n=3 was applied to each of these intervals.

Monitoring showed (see item 4) that L2- error for I1 is 1.×103, and for I2 is 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.185 and 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 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]).


  1. 1. A. Zygmund, Trigonometric Series, Vol. I, Cambridge, at the University Press, 1959
  2. 2. A. N. Krylov, About approximate calculations. Lectures given in 1906 (in Russian), St. Petersburg, Typolithography of K. Birkenfeld, 1907
  3. 3. A. N. Krylov, Lectures on approximate computing (in Russian), St. Petersburg, Printing house of Yu.N. Erlikh, 1911
  4. 4. K. S. Eckhoff, Accurate and efficient reconstruction of discontinuous functions from truncated series expansions, Math. Comp. v. 61,N204, 745–763, 1993
  5. 5. K. S. Eckhoff, Accurate reconstruction of functions of finite regularity from truncated Fourier series expansions, Math. Comp. v. 64,N210, 671–690, 1995
  6. 6. C. Lanczos, Discourse of Fourier Series, Edinburgh, Oliver and Boyd, 1966
  7. 7. W. B. Jones and G. Hardy, Accelerating convergence of trigonometric approximations, Math. Comp., 24, 547–560, 1970
  8. 8. A. Nersessian and A. Poghosyan, Accelerating the Convergence of Trigonometric Series, Central European Journal of Math., vol. 4:3,435–448, 2006
  9. 9. A. Barkhudaryan, R. Barkhudaryan and A. Poghosyan, Asymptotic Behavior of Eckhoff’s Method for Fourier Series Convergence Acceleration, Analysis in Theory and Applications, vol. 23 (3), 1–15, 2007
  10. 10. B. Adcock, Gibbs phenomenon and its removal for a class of orthogonal expansions, BIT, vol. 51(1), 7–41, 2011
  11. 11. B. Adcock, Gibbs phenomenon and its removal for a class of orthogonal expansions, BIT, vol. 51(1), 7–41, 2011
  12. 12. A. Poghosyan, On an auto-correction phenomenon of the Krylov-Gottlieb-Eckhoff method, IMA Journ. of Numer. Analysis, v. 31 (2), pp. 512–527, 2011
  13. 13. A. Nersessian, Quasi-polynomials of Bernoulli type and acceleration of convergence of Fourier series (in Russian). Reports of the National Academy of Sciences of Armenia, vol.104 (5), pp.280–286, 2004
  14. 14. A. Nersessian, Acceleration of convergence of Fourier-Bessel series for piecewise smooth functions (in Russian)
  15. 15. A. Nersessian, Accelerated convergence of eigenfunction expansions (in Russian). Reports of the National Academy of Sciences of Armenia, vol 107 (2), pp. 124–131,2007
  16. 16. Anry Nersessian, On an Over-Convergence Phenomenon for Fourier series. Basic Approach. Armen. J. Math., V. 10, N. 9, pp. 1–22, 2018
  17. 17. Anry Nersessian, A correction to the article “On an Over-Convergence Phenomenon for Fourier series. Basic Approach”. Armen. J. Math., V.11, N. 2 , pp. 1–2, 2019
  18. 18. A. Nersessian, Fourier tools are much more powerful than commonly thought. Lobachevskii Journal of Mathematics, vol.40, N.8 , pp.1122–1131, 2019
  19. 19. S. Wolfram, The Mathematica Book, fifth edition, Wolfram Media, 2003


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

Written By

Anry Nersessian

Submitted: 16 June 2020 Reviewed: 24 September 2020 Published: 02 December 2020