A Fast Method for Numerical Realization of Fourier Tools

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.


Introduction
One of the classical tools of Mathematics is the apparatus of Fourier series, based on the orthogonal system e iπkx ÈÉ , k ¼ 0, AE 1, AE 2, … complete in L 2 À1, 1 ½ . 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 À1, 1 ½ ), 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 À1, 1 ½ . 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 isgivenonthesegment À1, 1 ½ with Fourier coefficients f s ÈÉ , s ¼ 0, AE 1, … , AE n, n ≥ 1, andwiththefollowingjumpsofthefunctionf or its derivatives until degrees of the order q ≥ 1 at the points a k fg , À1 < a 1 < ⋯ < a m ¼ 1, 1 ≤ m < ∞, A p,k f ðÞ ¼ f k ðÞ a p À 0 ÀÁ À f k ðÞ a p þ 0 ÀÁ , In the neighborhoods of other points we assume that f ∈ C qþ1 . Let us construct a function g ¼ gx ðÞ , x ∈ À1, 1 ½ with Fourier coefficients g s ÈÉ ,0≤|s|≤n, which has the same jumps at the same points, and g ∈ 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 ¼ f À g ðÞ is q times continuously differentiable on whole axis, and then fx ðÞ À gx ðÞ¼ where r s x ðÞ¼os Àq ðÞ , s ! ∞, x ∈ À1, 1 ½ . Therefore, taking into account only first 2n þ 1 ðÞ Fourier coefficients and truncating the remainder term r s , it is possible to approximate f in the form However, the computing of the jumps A sk f ðÞ fg 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 f s ÈÉ (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 ½ 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 fg are given explicitly, converges to f with the rate on Àq ðÞ , n ! ∞.
K. Eckhoff suggests to find approximate values of jumpsÃ p,k ≈A p,k ÈÉ 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, … , mqþ 1 ðÞ , θn ≤|s k |≤n, 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].

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
where are Fourier coefficients of f , and D n ¼ d k fg , k ¼ 1, … , n, is a set of n different integers (n ≥ 1). We will assume that D 0 ¼ Ø.
Definition 2. Let n ≥ 1 be a fixed integer. Consider a system of functions where λ k fg are arbitrary parameters. Consider the linear span Q n ¼ span U n fg . We call a function q ∈ Q n as quasipolynomial of degree at most n.
It is easy to see that q ∈ Q n , if and only if when either qx ðÞ0orqx ðÞ¼ The number m will be considered below as the degree of the quasi-polynomial q. Definition 3. We call the system λ k fg ⊂  as parameters of the quasi-polynomial q x ðÞ¼ P k P β k x ðÞ exp i πλ k x ðÞ ∈ Q n and the number 1 þ β k as multiplicity of the parameter λ k . Definition 4. Let parameters λ k fg , k ¼ 1, … , n are 1 fixed. We denote by Q n λ k fg ðÞ the set of all corresponding q x ðÞ ∈ Q n . Remark 1. The set of quasi-polynomials q x ðÞ ∈ Q n is invariant with respect to a linear change of the variable x. The set of quasi-polynomials q x ðÞ ∈ Q n λ k fg ðÞ is invariant with respect to a shift of the variable x.
For a fixed integer n ≥ 1, consider a set of non-integer parameters λ k fg ⊂ , k ∈ D n , and the following infinite sequence (3) is invariant with respect to the numbering of the parameters λ k fg . Here this numbering is tied to the set k ∈ D n . In the general case, the parameters λ k fg (i.e., and t r,s ) may depend on n. In order not to complicate the notation, we will indicate this dependence as needed.

Remark 2. The last product in
Further we denote (see [17] for details) The system T r x ðÞ ,1=2exp i π rx ðÞ fg , r ∈ D n , is biorthogonal on the segment x ∈ À1, 1 ½ and L 2 -error of approximation fx ðÞ ≃ F n x ðÞ 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 (x ∈ À1, 1 ½ ) where sinc z ðÞ¼sin z ðÞ =z, sinc 0 ðÞ ¼ 1, z ∈ , plays a key role in the future. This Fourier series can be repeatedly differentiated by the parameter λ.

The case of integer or multiple parameters
To include the case of a quasi-polynomial q ∈ Q n containing some integer parameters from λ k fg , consider the following quasi-polynomials associated with formula (12) and sequence (9) Lemma 1. If λ j is an integer parameter, then where B k x ðÞis a Bernoulli polynomial (see (4)). Since the Bernoulli polynomial B k x ðÞ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
, k ≥ 1 can be decomposed into a Fourier series with coefficients (8), for which If in the sequence (9) parameters λ k fg and its corresponding multiplicities n k fg are known, then we will use the representation s ¼ 0, AE1, … ðÞ where r ∈ D n , D m ⊂ D n , n q ÈÉ are corresponding positive integers, and P j ∈ D m n j ¼ n, λ p 6 ¼ λ q if p 6 ¼ q. However, this sequence is still defined only for noninteger 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 j ∉ D n , λ j is integer, then t r,λ j ¼ ∞. So here the system T r fg from (10) does not exist. If λ j are integers for j, j ∈ D n , then, firstly, it is natural to accept t r,r ¼ 1, and secondly, we notice that for s 6 ¼ 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 ∈ . Bearing in mind Lemma 1, it is not difficult to proof, that. Lemma 2. Let in (9 * ) the parameters λ j ÈÉ , j ∈ D m , be given and the subset Λ consists of all integer parameters from λ j ÈÉ . Then.
i. If there is an λ k ∈ Λ, λ k ∉ D n , then the system T r fg from (10) does not exist. ii.

contains only different integers and Λ 2 ⊂ D m contains only multiple integers, then
T r x ðÞ¼exp i π rx ðÞ , r ∈ Λ 1 : With a λ r ∈ Λ 2 we can apply Lemma 1 and also find the explicit form of functions T r x ðÞ fg .
iii. If Λ 6 ¼ Ø and Λ 6 ¼ D n then the system T r x ðÞ fg , r ∈ D m nΛ is determined based on the sequence, similar to (15), in which n j ¼ 1, ∀j, and n is replaced by n À p ðÞ where p is the number of elements Λ.

Explicit form of the system T r fg
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 fg . Theorem 1. Suppose the sequence (15) is given and the possible integer parameters in λ k fg satisfy the condition λ k ∈ D m . Then the corresponding functions T r fg are quasi-polynomials and have the following explicit form where c r,j,k ¼ À1

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 A.
Let S n x ðÞbe a truncated Fourier series (see (7)). Consider formulas (9)-(10) with symbolic (not numerical) parameters λ q ÈÉ , q ∈ D n . Let us now choose a new set of integer numbersD n ¼d k To determine the parameters λ q ÈÉ , q ∈ D n , we additionally use Fourier coefficients f s ÈÉ , s ∈D 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 T r fg from (4) is used to approximate f by F n (see (10)). Here we show how Algorithm A can be step-by-step realized using given 2n ⋃ fd k no no (see (7)).
Step 1. Using the representation (9), we formally bring the left side of (18) to a common denominator, and fix the conditions λ q 6 ¼d k , k ¼ 1, 2, … , n. Then Eq. (18) will take the form (s ∈D n ) Step 2. Using the Vieta's formula decompose the products in (18) containing the parameters λ k fg 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 e k fg .
Step 3. Solve resulting equation by the least square method.
Step 4. Again, according to the Vieta's formulae find all roots z k fg , k ¼ 1, 2, … , n of the corresponding polynomial Pz ðÞ¼ P n p¼0 e p z nÀp , e 0 ¼ 1, and put λ k fg ¼ z k fg .

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 ∈ Q n (see Definition 2) and the sets D n ,D n and the Fourier coefficients f s ÈÉ , s ∈ D n ∪D n , of the f be given. Denote by Λ the set of integer parameters in the representation f ¼ P k P m k x ðÞexp i πμ k x ðÞ . In order for the approximation by Algorithm A to be exact (that is f x ðÞF n x ðÞ , x ∈ À1, 1 ½ ), it is necessary and sufficient that Λ ⊆ D n ∪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 ∈ L 2 À1, 1 ½ and the sets D n ,D n and the Fourier coefficients f s ÈÉ , s ∈ D n ∪D n , of the f be given. Applying Algorithm 1, we get parameters λ k fg , k ∈ D n . The accelerating the convergence of truncated Fourier series we define by formula (17).

Remark 5.
There are no Fourier coefficients f r ÈÉ , r ∈D m in formula (17). These coefficients are used in Algorithm A to an optimal choice of parameters λ k fg , k ∈ D m .

An numerical algorithm for Fourier transform
Consider the Fourier transform F of a function f ∈ L 2 À∞, ∞ ðÞ the form

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 On the other hand, our method allows one to approximate function f ≃ P f r T r on each segment k À h, k þ 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 T r x ðÞis a linear combination of functions of the form exp iπλ 1 x ðÞ or its derivative with respect to λ of a certain order. For example, if k ¼ 0, then for functions exp i πλ 1 x ðÞ , we have (see (12), here k ¼ 0) On the complex plane s ∈ , 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 in formula (21), interval À∞, ∞ ðÞ 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 .
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 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.

A simple example
Let us explain the details on the following example of an approximate finding of the Fourier transform (see notation (21)). and when deleting these semi-infinite intervals, further actions should lead to an error no greater than 5:99 Â 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 ¼À 2, 1=3 ðÞ , I 2 ¼ 1=3, 3=2 ðÞ . 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: Â 10 À3 , and for I 2 is 1:44 Â 10 À3 . We see that this is less than allowed above.
The total L 2 -error turned out to be 4:38 Â 10 À3 . As for the uniform error, it is reached approximately at point ω ¼ 0:185 and is equal to 1:04 Â 10 À3 .
These calculations were performed using system Wolfram Mathematica 9.

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]).