A rational approximation of the Fourier transform by integration with exponential decay multiplier

Recently we have reported a new method of rational approximation of the sinc function obtained by sampling and the Fourier transforms. However, this method requires a trigonometric multiplier that originates from shifting property of the Fourier transform. In this work we show how to represent the Fourier transform of a function $f(t)$ in form of a ratio of two polynomials without any trigonometric multiplier. A MATLAB code showing algorithmic implementation of the proposed method for rational approximation of the Fourier transform is presented.


Introduction
The forward and inverse Fourier transforms of two related functions f (t) and F (ν) can be defined in a symmetric form as [1,2] and where variables t and ν are the corresponding Fourier-transformed arguments in t-space and ν-space, respectively (time t vs. frequency ν, for example). Fourier transform methods are widely used in many applications including signal processing [1,2], spectroscopy [3,4] and computational finance [5][6][7].
There are several efficient methods have been reported for rational approximations in literature. For example, the rational approximations may be built on the basis of the Newman nodes [8], Chebyshev nodes [9], logarithmic nodes [10] and so on.
Recently we have reported a new method of rational approximation of the Fourier transform (1) as given by [11] F {f (t)} (ν) ≈ e 2πiνa where M is an integer determining number of summation terms 2 M −1 , a is a shift constant, σ is a decay (damping) constant and It has been noticed that approximation (3) is not purely rational and there was a question whether or not a rational function of the Fourier transform (1) in explicit form without any trigonometric multiplier of kind e 2πiνa = cos (2πνa) + i sin (2πνa) , depending on argument ν, can be obtained [12]. Theoretical analysis shows that this trigonometric multiplier originating from shifting property of the Fourier transform can be indeed excluded. As a further development of our work [11], in this paper we derive a rational function of the Fourier transform (1) that has no any trigonometric multiplier of kind (4). Therefore, it can be used as an alternative to the Padé approximation. To the best of our knowledge, this method of rational approximation of the Fourier transform (1) for a non-periodic function f (t) has never been reported in scientific literature.

Preliminaries
Assume that Re {f (t)} is even while Im {f (t)} is odd such that f : R → C, but Re {f } : R → R and Im {f } : R → R. Then it is not difficult to see that the Fourier transform (1) of the function f (t) can be expanded into two integral terms as follows Assume also that the function f (t) behaves in such a way that for some positive numbers τ 1 and τ 2 the following integrals are negligibly small and can be ignored in computation. Consequently, we can approximate the Fourier transform as given by Further the values 2τ 1 and 2τ 2 will be regarded as widths (pulse widths) for the real and imaginary parts of the function f (t), respectively.

New sampling method
Consider a sampling formula (see, for example, equation (3) in [13]) where is the sinc function, t n is a set of sampling points, h is small adjustable parameter (step) and ε (t) is error term. François Viète discovered that the sinc function can be represented by cosine product 1 [14,15] In our earlier publications we introduced a product-to-sum identity [16] M m=1 cos and applied it for sampling [17,18] as incomplete cosine expansion of the sinc function for efficient computation of the Voigt/complex error function. It is worth noting that this product-to-sum identity has also found some efficient applications in computational finance [6,19,20] involving numerical integration.
Change of variable 2 M −1 → M in the limit above leads to Therefore, by truncating integer M and by making another change of variable t → πt/h we obtain The right side of equation (9) is periodic due to finite number of the summation terms. As a result, the approximation (9) is valid only within the At equidistantly separated sampling grid-points such that t n = nh, the substitution of approximation (9) into sampling formula (6) gives It is important that in sampling procedure the total number of the sampling grid-points 2N + 1 as well as the step h should be properly chosen to insure that the widths 2τ 1 and 2τ 2 are entirely covered. As we can see, the sampling formula (10) is based on incomplete cosine expansion of the sinc function that was proposed in our previous works [17,18] as a new approach for rapid and highly accurate computation of the Voigt/complex error function [21][22][23]. Computations we performed show that this method of sampling is particularly efficient in numerical integration.

Even function
Suppose that our objective is to approximate the sinc function sinc (πν). First we take the inverse Fourier transform (2) of the sinc function is known as the rectangular function. This function is even since rect (t) = rect (−t). The rectangular function rect (t) has two discontinuities at t = −1/2 and t = 1/2. Therefore, it is somehow problematic to perform sampling over this function. However, we can use the fact that Thus, by taking a sufficiently large value for the integer k, say k = 35, we can approximate the rectangular function (11) quite accurately as Figure 1 shows the function f (t) = 1/ (2t) 70 + 1 by blue curve. As we can see from this figure, the function very rapidly decreases at |t| > 1/2 with increasing t. Therefore, we can take τ 1 = 0.6. Thus, the width of this function is 2τ 1 = 1.2.
Sampling of function f (t) = 1/ (2t) 70 + 1 in accordance with equation (10) results in a periodic dependence. Consequently, due to periodicity on the right side of equation (10) it cannot be utilized for rational approximation of the Fourier transform. However, this problem can be effectively resolved by sampling the function f (t) e σt instead of f (t) itself. This leads to   (12) in form can effectively eliminate this periodicity due to presence of the exponential decay multiplier e −σt on the right side. This suppression effect can be seen from the Fig. 3 illustrating the results of computation for the even function f (t) = 1/ (2t) 70 + 1 by approximation (13) at M = 32, N = 28 with σ = 0 (blue curve), σ = 0.25 (red curve) and σ = 0.75 (green curve). As it is depicted by blue curve, at σ = 0 the function is periodic. However, as decay coefficient σ increases, the exponential multiplier e −σt suppresses all the peaks (except the first peak at the origin) such that the resultant function tends to become solitary along the entire positive t-axis. This tendency can be observed by red and green curves at σ = 0.25 and σ = 0.75, respectively. As a consequence, if the damping multiplier σ is big enough, say greater than Thus, substituting approximation (13) into equation (5) and considering the fact that at sufficiently large σ the function becomes solitary along positive x-axis, the upper limit τ 1 of integration can be replaced by infinity This integral can be taken analytically in form of rational function now and after some trivial rearrangements that exclude double summation, it follows where the expansion coefficients are given by

Odd function
Consider, as an example, the following function We can see that the condition t rect (t) = − (−t rect (−t)) for odd function in its imaginary part is satisfied. The function Im {f (t)} = t/ (2t) 70 + 1 is shown in the Fig. 1 by red curve. We can take τ 2 = 0.6 and the width is 2τ 2 = 1.2.
Using exactly same procedure as it has been described above and considering the fact that at sufficiently large σ the upper limit τ 2 of integration can be replaced by infinity, we can write 3

This leads to
where the expansion coefficients are The Fourier transform of the function it rect (t) can be readily found analytically Gray curve in Fig. 4 illustrates the Fourier transform of the function f (t) = it/ (2t) 70 + 1 obtained by using approximation (15) at M = 32, N = 28, h = 0.04 and σ = 3. The original function is also shown for comparison by black dashed curve. These two curves in the interval −2π ν 2π are also not distinctive visually.
In our recent work [24] we applied alternative method of sampling by using incomplete cosine expansion of the Gaussian function of kind h e −(t/c) 2 / (c √ π), where c and h are the fitting parameters. We have shown that this method of sampling can also be used to obtain high-accuracy computation of the Voigt/complex error function. In our future work will apply this method of sampling as an alternative that may reduce the absolute difference for rational approximations of the piecewise functions with discontinuities.

Alternative representation
For a function f (t) = Re {f (t)} + i Im {f (t)}, where its real part Re {f (t)} is even and its imaginary part Im {f (t)} is odd, we can write Using a Computer Algebra System (CAS) supporting symbolic programming it is not difficult to find coefficients p k and q k to represent this approximation as where are polynomials of the orders 4M − 1 and 4M , respectively. Padé approximation is one of the efficient methods to represent a function in form of ratio of two polynomials. Our preliminary numerical results show that the proposed new method of rational approximation may significantly extend the range [ν min , ν max ] in coverage [12] than the conventional Padé approximation.

MATLAB code and description
The MATLAB code shown below is written as a function file raft.m that can be simply copied and pasted to create m-file in the MATLAB environment. The name of this function file originates from the abbreviation RAFT that stands for Rational Approximation of the Fourier Transform. The command raft(opt) performs sampling and then computation of the expansion coefficients α m , β m , η m , θ m , κ m , λ m , µ m . Once the coefficients are determined, the program executes the Fourier transform according to equations (14) and (15) This method of the rational approximation is based on integration involving exponential decay multiplier e −σt . The computational test we performed shows that this method of the Fourier transform can provide relatively accurate approximations of the Fourier transform even for the functions with discontinuities like rect (t) and it rect (t). Furthermore, this method shows that for the well-behaved function f (t) = √ πe −(πt) 2 + i π 3/2 te −(πt) 2 with only 16 summation terms the rational approximations (14) and (15) provide the Fourier transform with absolute differences not exceeding 3 × 10 −10 and 9 × 10 −10 for its real and imaginary parts, respectively. Our preliminary results indicate that the proposed method may be promising for rational approximation over the wide range [ν min , ν max ].