A rational approximation of the Fourier transform by integration with exponential decay multiplier
AA rational approximation of the
Fourier transform by integration withexponential decay multiplier
S. M. Abrarov
1, 2 , R. Siddiqui
1, 2, 3 , R. K. Jagpal
2, 3 , and B. M. Quine
1, 31Dept. Earth and Space Science and Engineering, York University, 4700 Keele St., Canada, M3J 1P32Epic College of Technology, 5670 McAdam Rd., Mississauga, Canada, L4Z 1T23Dept. Physics and Astronomy, York University, 4700 Keele St., Toronto, Canada, M3J 1P3
January 15, 2020
Abstract
Recently we have reported a new method of rational approximationof the sinc function obtained by sampling and the Fourier transforms.However, this method involves a trigonometric multiplier e πiνa = cos (2 πνa ) + i sin (2 πνa ) , where ν is the argument and a is the shift constant. In this work weshow how to avoid this trigonometric multiplier in order to representthe Fourier transform of function f ( t ) in explicit form as a rationalapproximation of kind F { f ( t ) } ( ν ) = (cid:90) ∞−∞ f ( t ) e − πiνt dν ≈ P ( ν ) Q ( ν ) , where P ( ν ) and Q ( ν ) are polynomials. A MATLAB code showingalgorithmic implementation of the proposed method for rational ap-proximation of the Fourier transform is presented. Keywords: rational approximation; Fourier transform; sampling;sinc function a r X i v : . [ m a t h . G M ] J a n Introduction
The forward and inverse Fourier transforms of two related functions f ( t ) and F ( ν ) can be defined in a symmetric form as [7, 10] F { f ( t ) } ( ν ) = F ( ν ) = ∞ (cid:90) −∞ f ( t ) e − πiνt dt (1)and F − { F ( ν ) } ( t ) = f ( t ) = ∞ (cid:90) −∞ F ( ν ) e πiνt dν, (2)where variables t and ν are the corresponding Fourier-transformed argumentsin t -space and ν -space, respectively (time t vs. frequency ν , for example).Recently we have reported a new method of rational approximation of theFourier transform (1) as given by [4] F { f ( t ) } ( ν ) ≈ e πiνa M − (cid:88) m =1 A m ( σ + 2 πiν ) + B m C m + ( σ + 2 πiν ) , (3)where M is an integer determining number of summation terms 2 M − , a is ashift constant, σ is a decay (damping) constant and A m = 12 M − N (cid:88) n =0 f ( nh − a ) e σnh cos ( C m nh ) ,B m = 12 M − N (cid:88) n =0 f ( nh − a ) e σnh C m sin ( C m nh ) ,C m = π (2 m − M h . are expansion coefficients.It has been noticed that approximation (3) is not purely rational and therewas a question whether or not a rational function of the Fourier transform(1) in explicit form without any trigonometric multiplier of kind e πiνa = cos (2 πνa ) + i sin (2 πνa ) , (4)2epending on argument ν , can be obtained [15]. Theoretical analysis showsthat this trigonometric multiplier can be indeed excluded. In this work wederive a rational function of the Fourier transform (1) that has no any trigono-metric multiplier of kind (4). The computational test we performed revealsthat this method of rational approximation of the Fourier transform may bemore accurate in computation as compared to that of based on approxima-tion (3). To the best of our knowledge, this method of rational approximationof the Fourier transform (1) has never been reported in scientific literature. 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 seethat the Fourier transform (1) of the function f ( t ) can be expanded into twointegral terms as follows F { f ( t ) } ( ν ) = F ( ν ) = 2 ∞ (cid:90) Re { f ( t ) } cos (2 πνt ) dt + 2 ∞ (cid:90) Im { f ( t ) } sin (2 πνt ) dt. Assume also that the function f ( t ) behaves in such a way that for somepositive numbers τ and τ the following integrals ∞ (cid:90) τ Re { f ( t ) } cos (2 πνt ) dt ≈ ∞ (cid:90) τ Im { f ( t ) } sin (2 πνt ) dt ≈ F { f ( t ) } ( ν ) = F ( ν ) ≈ τ (cid:90) Re { f ( t ) } cos (2 πνt ) dt + 2 τ (cid:90) Im { f ( t ) } sin (2 πνt ) dt. (5) Further the values 2 τ and 2 τ will be regarded as widths (pulse widths) forthe real and imaginary parts of the function f ( t ), respectively.3 .2 New sampling method Consider a sampling formula (see, for example, equation (3) in [14]) f ( t ) = N (cid:88) n = − N f ( t n ) sinc (cid:16) πh ( t − t n ) (cid:17) + ε ( t ) , (6)where sinc ( t ) = sin tt , t (cid:54) = 01 , t = 0 , is the sinc function, t n is a set of sampling points, h is small adjustableparameter (step) and ε ( t ) is error term. Fran¸cois Vi`ete discovered that thesinc function can be represented by cosine product [9, 11]sinc ( t ) = ∞ (cid:89) m =1 cos (cid:18) t m (cid:19) . (7)In our earlier publications we introduced a product-to-sum identity [16] M (cid:89) m =1 cos (cid:18) t m (cid:19) = 12 M − M − (cid:88) m =1 cos (cid:18) m − M t (cid:19) (8)and applied it for sampling [2, 3] as incomplete cosine expansion of the sincfunction for efficient computation of the Voigt/complex error function. Itis worth noting that this product-to-sum identity has also found some effi-cient applications in Computational Finance [8, 12, 13] involving numericalintegration.Comparing identities (7) and (8) immediately yieldssinc ( t ) = lim M →∞ M − M − (cid:88) m =1 cos (cid:18) m − / M − t (cid:19) . Unlike equation (7), this limit consists of sum of cosines instead of product ofcosines. As a result, its application provides significant flexibilities in variousnumerical integrations [2, 3, 8, 12, 13]. This equation is also attributed to Euler. M − → M in the limit above leads tosinc ( t ) = lim M →∞ M M (cid:88) m =1 cos (cid:18) m − / M t (cid:19) . Therefore, by truncating integer M and by making another change of variable t → πt/h we obtainsinc (cid:16) πh t (cid:17) ≈ M M (cid:88) m =1 cos (cid:18) π ( m − / M h t (cid:19) , − M h (cid:54) t (cid:54) M h. (9)The right side of equation (9) is periodic due to finite number of the sum-mation terms. As a result, the approximation (9) is valid only within theinterval t ∈ [ − M h, M h ] . At equidistantly separated sampling grid-points such that t n = nh , thesubstitution of approximation (9) into sampling formula (6) gives f ( t ) ≈ M M (cid:88) m =1 N (cid:88) n = − N f ( nh ) cos (cid:18) π ( m − / M h ( t − nh ) (cid:19) , − M h (cid:54) t (cid:54) M h. (10)
It is important that in sampling procedure the total number of the samplinggrid-points 2 N + 1 as well as the step h should be properly chosen to insurethat the widths 2 τ and 2 τ are entirely covered.As we can see, the sampling formula (10) is based on incomplete co-sine expansion of the sinc function that was proposed in our previous works[2, 3] as a new approach for rapid and highly accurate computation of theVoigt/complex error function [1, 5, 6]. Computations we performed show thatthis method of sampling is particularly efficient in numerical integration. Suppose that our objective is to approximate the sinc function sinc ( πν ).First we take the inverse Fourier transform (2) of the sinc function F − { sinc ( πν ) } ( t ) = ∞ (cid:90) −∞ sinc ( πν ) e − πiνt dν = rect ( t ) , where rect ( t ) = , if | t | < / / , if | t | = 1 / , if | t | > / ,
5s known as the rectangular function. This function is even since rect ( t ) =rect ( − t ). The rectangular function rect ( t ) has two discontinuities at t = − / t = 1 /
2. Therefore, it is somehow problematic to perform samplingover this function. However, we can use the fact thatrect ( t ) = lim k →∞ t ) k + 1 . (11)Thus, by taking a sufficiently large value for the integer k , say k = 35, wecan approximate the rectangular function (11) quite accurately asrect ( t ) ≈ f ( t ) = 1(2 t ) + 1 . - - t - /[( ) + ] , t /[( ) + ] Fig. 1.
The even / (cid:0) (2 t ) + 1 (cid:1) and odd t/ (cid:0) (2 t ) + 1 (cid:1) functionsshown by blue and red curves, respectively. Figure 1 shows the function f ( t ) = 1 / (cid:0) (2 t ) + 1 (cid:1) by blue curve. Aswe can see from this figure, the function very rapidly decreases at | t | > / t . Therefore, we can take τ = 0 .
6. Thus, the width of thisfunction is 2 τ = 1 . f ( t ) = 1 / (cid:0) (2 t ) + 1 (cid:1) in accordance with equation(10) results in a periodic dependence. Consequently, due to periodicity onthe right side of equation (10) it cannot be utilized for rational approximation6f the Fourier transform. However, this problem can be effectively resolvedby sampling the function f ( t ) e σt instead of f ( t ) itself. This leads to f ( t ) e σt ≈ M M (cid:88) m =1 N (cid:88) n = − N f ( nh ) e σnh cos (cid:18) π ( m − / M h ( t − nh ) (cid:19) , − M h (cid:54) t (cid:54) M h. (12) ν - - - Fig. 2.
Approximation (12) to the function f ( t ) e σt = e σt / (cid:2) (2 t ) + 1 (cid:3) computed at M = 32 , N = 28 , h = 0 . with σ = 0 (blue curve), σ = 0 . (red curve) and σ = 0 . (greencurve). Figure 2 shows the results of computation for even function f ( t ) =1 / (cid:0) (2 t ) + 1 (cid:1) by approximation (12) at M = 32, N = 28, h = 0 .
04 with σ = 0 (blue curve), σ = 0 .
25 (red curve) and σ = 0 .
75 (green curve). As wecan see from this figure, all three curves are periodic as expected. However,if the constant σ is big enough, then slight rearrangement of equation (12)in form f ( t ) ≈ e − σt M M (cid:88) m =1 N (cid:88) n = − N f ( nh ) e σnh cos (cid:18) π ( m − / M h ( t − nh ) (cid:19) , (13)can effectively eliminate this periodicity due to presence of the exponentialdecay multiplier e − σt on the right side. This suppression effect can be seenfrom the Fig. 3 illustrating the results of computation for the even function7 ( t ) = 1 / (cid:0) (2 t ) + 1 (cid:1) by approximation (13) at M = 32, N = 28 with σ = 0 (blue curve), σ = 0 .
25 (red curve) and σ = 0 .
75 (green curve). Asit is depicted by blue curve, at σ = 0 the function is periodic. However, asdecay coefficient σ increases, the exponential multiplier e − σt suppresses allthe peaks (except the first peak at the origin) such that the resultant functiontends to become solitary along the entire positive t -axis. This tendency canbe 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 thanunity, the approximated function becomes practically solitary as the originalfunction f ( t ) = 1 / (cid:0) (2 t ) + 1 (cid:1) itself.Thus, substituting approximation (13) into equation (5) and consider-ing the fact that at sufficiently large σ the function becomes solitary alongpositive x -axis, the upper limit τ of integration can be replaced by infinityas F { f ( t ) } ( ν ) = F (cid:40) t ) + 1 (cid:41) ( ν ) ≈ τ (cid:90) (cid:34) e − σt M M (cid:88) m =1 N (cid:88) n = − N f ( nh ) e σnh cos (cid:18) π ( m − / M h ( t − nh ) (cid:19)(cid:35) cos (2 πνt ) dt ≈ ∞ (cid:90) (cid:34) e − σt M M (cid:88) m =1 N (cid:88) n = − N f ( nh ) e σnh cos (cid:18) π ( m − / M h ( t − nh ) (cid:19)(cid:35) cos (2 πνt ) dt. This integral can be taken analytically in form of rational function now andafter some trivial rearrangements that exclude double summation, it followsthat
F { f ( t ) } ( ν ) ≈ M (cid:88) m =1 α m + β m ν κ m + λ m ν + ν , (14)where the expansion coefficients are given by α m = 18 M π N (cid:88) n = − N f ( nh ) e nhσ (cid:0) µ m + σ (cid:1) ( σ cos ( nhµ m ) + µ m sin ( nhµ m )) , β m = 12 M π N (cid:88) n = − N f ( nh ) e nhσ ( σ cos ( nhµ m ) − µ m sin ( nhµ m )) , For this integration we imply that the interval 2
N h along t -axis occupied by samplinggrid-points is larger than the function width 2 τ = 1 . ν - - Fig. 3.
Evolution to the function f ( t ) = 1 / (cid:2) (2 t ) + 1 (cid:3) computedby approximation (13) at M = 32 , N = 28 , h = 0 . with σ = 0 (blue curve), σ = 0 . (red curve) and σ = 0 . (green curve). κ m = 116 π (cid:0) µ m + σ (cid:1) ,λ m = 12 π (cid:0) σ − µ m (cid:1) and µ m = π ( m − / M h .
Figure 4 shows the original sinc function sinc ( ν ) and its approximation(14) within the interval − π (cid:54) ν (cid:54) π at M = 32 , N = 28, h = 0 .
04 and σ = 2 .
75 by black dashed and light blue curves, respectively. These twocurves are not visually distinctive.
Consider, as an example, the following function f ( t ) = it (2 t ) + 1 ≈ it rect ( t ) . We can see that the condition t rect (t) = − ( − t rect ( − t)) for odd functionin its imaginary part is satisfied. The function Im { f ( t ) } = t/ (cid:0) (2 t ) + 1 (cid:1) - - ν - Fig. 4.
Approximations of the functions sinc ( πν ) and (sin ( πν ) − πν cos ( πν )) / (cid:0) πν ) (cid:1) within interval − π ≤ ν ≤ π .Both approximations are obtained by equation (14) for input functions / (cid:0) (2 t ) + 1 (cid:1) and it/ (cid:0) (2 t ) + 1 (cid:1) at M = 32 , N = 28 , h = 0 . , σ = 2 . (light blue curve) and at M = 32 , N = 28 , h = 0 . , σ = 3 (gray curve), respectively. The original functions sinc ( πν ) and (sin ( πν ) − πν cos ( πν )) / (cid:0) πν ) (cid:1) are also shown by black dashedcurves for comparison. is shown in the Fig. 1 by red curve. We can take τ = 0 . τ = 1 . σ the upper limit τ of integration canbe replaced by infinity, we can write F { f ( t ) } ( ν ) = F (cid:40) it (2 t ) + 1 (cid:41) ( ν ) ≈ τ (cid:90) (cid:34) e − σt M M (cid:88) m =1 N (cid:88) n = − N f ( nh ) e σnh cos (cid:18) π ( m − / M h ( t − nh ) (cid:19)(cid:35) sin (2 πνt ) dt ≈ ∞ (cid:90) (cid:34) e − σt M M (cid:88) m =1 N (cid:88) n = − N f ( nh ) e σnh cos (cid:18) π ( m − / M h ( t − nh ) (cid:19)(cid:35) sin (2 πνt ) dt. In this integration we imply again that the interval 2
N h along t -axis occupied bysampling grid-points is larger than the function width 2 τ = 1 . F { f ( t ) } ( ν ) ≈ M (cid:88) m =1 η m ν + θ m ν κ m + λ m ν + ν , (15)where the expansion coefficients are η m = 14 M π N (cid:88) n = − N f ( nh ) e nhσ (cid:0)(cid:0) σ − µ m (cid:1) cos ( nhµ m ) + 2 σµ m sin ( nhµ m ) (cid:1) and θ m = 1 M π N (cid:88) n = − N f ( nh ) e nhσ cos ( nhµ m ) . The Fourier transform of the function it rect ( t ) can be readily foundanalytically F { it rect ( t ) } ( ν ) = ∞ (cid:90) −∞ it rect ( t ) e − πiνt dt = / (cid:90) − / it rect ( t ) e − πiνt dt = / (cid:90) − / it e − πiνt dt = sin ( πν ) − πν cos ( πν )2 ( πν ) . Gray curve in Fig. 4 illustrates the Fourier transform of the function f ( t ) = it/ (cid:0) (2 t ) + 1 (cid:1) obtained by using approximation (15) at M = 32, N = 28, h = 0 .
04 and σ = 3. The original function F { it rect ( t ) } ( ν ) = sin ( πν ) − πν cos ( πν )2 ( πν ) is also shown for comparison by black dashed curve. These two curves in theinterval − π (cid:54) ν (cid:54) π are also not distinctive visually. Figures 5 shows the absolute difference between original sinc function sinc ( πν )and its approximation (14) for input function f ( t ) = 1 / (cid:0) (2 t ) + 1 (cid:1) calcu-lated at M = 32, N = 28, h = 0 .
04 and σ = 2 .
7. As we can see, the absolute11 - - ν Fig. 5.
Absolute difference between the original sinc function sinc ( πν ) and its rational approximation (14) for input function f ( t ) = 1 / (cid:0) (2 t ) + 1 (cid:1) at M = 32 , N = 28 , h = 0 . and σ = 2 . . difference within the interval − π (cid:54) ν (cid:54) π does not exceed 2 . × − . Thisaccuracy is better than that of shown in our recent publication, where weused equation (3) for the sinc function approximation (see Fig. 6 in [4]).Figure 6 shows the absolute difference between original function givenby (sin ( πν ) − πν cos ( πν )) / (cid:0) πν ) (cid:1) and its approximation (15) for inputfunction f ( t ) = it/ (cid:0) (2 t ) + 1 (cid:1) calculated at M = 32, N = 28, h = 0 . σ = 3. We can see that the absolute difference within the interval − π (cid:54) ν (cid:54) π does not exceed 6 × − .It should be noted that with more well-behaved functions we can obtainconsiderably higher accuracies. For example, suppose that we need to obtainthe Fourier transform of the function f ( t ) = √ πe − ( πt ) + i (cid:104) π / te − ( πt ) (cid:105) byusing approximations (14) and (15). Analytically, its Fourier transform is F (cid:110) √ πe − ( πt ) + i (cid:104) π / te − ( πt ) (cid:105)(cid:111) ( ν ) = e − ν + νe − ν , where e − ν is the Fourier transform of √ πe − ( πt ) while νe − ν is the Fouriertransform of iπ / te − ( πt ) .Blue curve in Fig. 7 corresponds to the absolute difference between func-tion e − ν and its approximation (14) for input function √ πe − ( πt ) at M = 16,12 - - ν Fig. 6.
Absolute difference between the original function (sin ( πν ) − πν cos ( πν )) / (cid:0) πν ) (cid:1) and its rational approximation (15) for input function f ( t ) = it/ (cid:0) (2 t ) + 1 (cid:1) at M = 32 , N = 28 , h = 0 . and σ = 3 . N = 23, h = 0 .
119 and σ = 6 .
9. Red curve in Fig. 7 corresponds to theabsolute difference between function νe − ν and its approximation (15) forinput function iπ / te − ( πt ) at M = 16, N = 23, h = 0 .
119 and σ = 5 .
9. Wecan see that with only 16 summation terms the absolute differences do notexceed 3 × − and 9 × − . These results demonstrate that the rationalapproximations (14) and (15) can be highly accurate in the Fourier transformof well-behaved functions. For a function f ( t ) consisting of both, the real and imaginary parts, we canwrite F { f ( t ) } ( ν ) ≈ M (cid:88) m =1 α m + β m ν κ m + λ m ν + ν + M (cid:88) m =1 η m ν + θ m ν κ m + λ m ν + ν or F { f ( t ) } ( ν ) ≈ M (cid:88) m =1 α m + η m ν + β m ν + θ m ν κ m + λ m ν + ν . - ν × - × - × - × - Fig. 7.
Absolute difference between the original functions e − ν and its approximation (14) for input function √ πe − ( πt ) at M = 16 , N = 23 , h = 0 . and σ = 6 . (blue curve). Absolute differencebetween the original functions νe − ν and its approximation (15) forinput function iπ / te − ( πt ) at M = 16 , N = 23 , h = 0 . and σ = 5 . (red curve). Using a Computer Algebra System (CAS) supporting symbolic programmingit is not difficult to find coefficients p k and q k to represent this approximationas F { f ( t ) } ( ν ) ≈ P ( ν ) Q ( ν ) , where P ( ν ) = p + p ν + p ν + · · · + p M − ν M − + p M − ν M − and Q ( ν ) = q + q ν + q ν + · · · + q M − ν M − + q M ν M , are polynomials of the orders 4 M − M , respectively. The MATLAB code shown below is written as a function file raft.m that canbe simply copied and pasted to create m -file in the MATLAB environment.14he name of this function file originates from the abbreviation RAFT thatstands for Rational Approximation of the Fourier Transform. The command raft(opt) performs sampling and then computation of the expansion coeffi-cients α m , β m , η m , θ m , κ m , λ m , µ m . Once the coefficients are determined, theprogram executes the Fourier transform according to equations (14) and (15)for even and odd input functions, respectively. The results of computationsare generated in two plots. The first plot shows the Fourier transform ofinput function while the second plot illustrates its absolute difference.There are four option values for opt argument. At opt = 0 , opt = 1 , opt= 2 and opt = 3 the corresponding input functions are rect (t), it rect (t), √ πe − ( πt ) and iπ / te − ( πt ) . The default value is opt = 0 signifying that forthe commands without argument raft and raft() , the value zero for opt isassigned.The authors did not attempt to optimize the code but rather to write itin a simple way with required comment lines in order to make it clear andintuitive for reading. The program was built and tested on MATLAB 2014a.However, the code should run in any version of MATLAB since it utilizes themost common commands. function raft(opt)% This function file performs the Rational Approximation of the Fourier% Transform (RAFT) for some functions and generates figures showing their% Fourier transform (FT) and absolute differencies.%% SYNOPSIS:% opt = 0 provides FT for input function: f(t) = rect(t)% opt = 1 provides FT for input function: f(t) = i*t*rect(t)% opt = 2 provides FT for input function:% f(t) = i*sqrt(pi)*exp(-(pi*t)^2)% opt = 3 provides FT for input function:% f(t) = i*pi^(3/2)*t*exp(-(pi*t)^2)%% The code is written by authors of this paper, York University, Toronto,% Canada, December 2019.clcif nargin == 0disp( ' Missing input parameter! Option opt = 0 is assigned. ' )opt = 0; ndswitch optcase 0M = 32; N = 28; h = 0.04; sigma = 2.7; % define parametersn = -N:N; % define array nf = 1./((2*n*h).^70 + 1); % 1st even functioncase 1M = 32; N = 28; h = 0.04; sigma = 3; % define parametersn = -N:N; % define array nf = n*h./((2*n*h).^70 + 1); % odd function, imaginary unit i is ...% omittedcase 2M = 16; N = 23; h = 0.119; sigma = 6.9; % define parametersn = -N:N; % define array nf = sqrt(pi)*exp(-(pi*n*h).^2); % 2nd even functioncase 3M = 16; N = 23; h = 0.119; sigma = 5.9; % define parametersn = -N:N; % define array nf = pi^(3/2)*n*h.*exp(-(pi*n*h).^2); % odd function, imaginary ...% unit i is omittedotherwisedisp( ' Wrong parameter! Enter either 0, 1, 2 or 3. ' )returnend% Initiate arrays for the expansion coefficientsalpha = zeros(M,1);beta = zeros(M,1);eta = zeros(M,1);theta = zeros(M,1);kappa = zeros(M,1);lambda = zeros(M,1);% Compute the expansion coefficientsf = f.*exp(n*h*sigma); % redefine the function arrayfor m = 1:Mmu = pi*(m - 1/2)/M; % must be defined firstif mod(opt,2) == 0 % if function is even, then:alpha(m) = 1/(8*M*pi^4)*sum(f.*((mu/h)^2 + sigma^2).*(sigma* ...cos(n*mu) + mu/(h)*sin(n*mu)));beta(m) = 1/(2*M*pi^2)*sum(f.*(sigma*cos(n*mu) - mu/h* ...sin(n*mu))); lse % otherwise, if function is odd then:eta(m) = 1/(4*M*pi^3)*sum(f.*((sigma^2 - (mu/(h))^2)* ...cos(n*mu) + 2*sigma*mu/h*sin(n*mu)));theta(m) = 1/(M*pi)*sum(f.*cos(n*mu));endkappa(m) = 1/(16*pi^4)*((mu/h)^2 + sigma^2)^2;lambda(m) = 1/(2*pi^2)*(sigma^2 - (mu/h)^2);end% Rational approximation (14) in string formatstrEven = [ ' f + (alpha(m) + beta(m)*nu.^2)./(kappa(m) + ' , ... ' lambda(m)*nu.^2 + nu.^4) ' ]; % for even function% Rational approximation (15) in string formatstrOdd = [ ' f + (eta(m)*nu + theta(m)*nu.^3)./(kappa(m) + ' , ... ' lambda(m)*nu.^2 + nu.^4) ' ]; % for odd functionnu = linspace(-2*pi,2*pi,1000); % define array for the argument nuf = 0; % reset function to zeroswitch optcase {0, 2} % two even functionsfor m = 1:M % main computation for even functionsf = eval(strEven); % use approximation (14)endif opt == 0 % assign some strings to plotdisp( ' Input function: f(t) = rect(t) ' )fStr = ' sinc(nu) ' ;yStr = ' \it{sinc(\pi\nu)} ' ;else % if opt = 2, then:disp( ' Input function: f(t) = sqrt(pi)*exp(-(pi*t)^2) ' )fStr = ' exp(-nu.^2) ' ;yStr = ' \it{exp(-\nu^2)} ' ;endcase {1, 3} % two odd functionsfor m = 1:M % main computation for odd functionsf = eval(strOdd); % use approximation (15)endif opt == 1 % assign some strings to plotdisp( ' Input function: f(t) = i*t*rect(t) ' )fStr = ' (sin(pi*nu) - pi*nu.*cos(pi*nu))./(2*(pi*nu).^2) ' ;yStr = ' \it{(sin(\pi\nu)-\pi\nu cos(\pi\nu))/(2(\pi\nu)^2)} ' ;else % if opt = 3, then:disp( ' Input function: f(t) = i*pi^(3/2)*t*exp(-(pi*t)^2} ' ) Str = ' nu.*exp(-nu.^2) ' ;yStr = ' \it{\nu exp(-\nu^2)} ' ;endendfunc2plot(nu,f,fStr,yStr); % call function to plotabsDiff(nu,f,fStr); % call to plot the abdolute differencefunction func2plot(nu,data2plot,funcStr,yAxisStr)% FIGURE 1figure1 = figure;axes1 = axes( ' Parent ' ,figure1, ' FontSize ' ,12);xlim(axes1,[-2*pi,2*pi]);box(axes1, ' on ' );grid(axes1, ' on ' );hold(axes1, ' all ' );plot1 = plot(nu,data2plot,nu,eval(funcStr), ' Parent ' ,axes1);set(plot1(1), ' LineWidth ' ,3, ' Color ' ,[0 1 1]);set(plot1(2), ' LineStyle ' , ' -- ' , ' LineWidth ' ,2);xlabel( ' Parameter \it{\nu} ' , ' FontSize ' ,14);ylabel(yAxisStr, ' FontSize ' ,14);endfunction absDiff(nu,data2plot,funcStr)% FIGURE 2figure2 = figure;axes2 = axes( ' Parent ' ,figure2, ' FontSize ' ,12);xlim(axes2,[-2*pi,2*pi]);box(axes2, ' on ' );grid(axes2, ' on ' );hold(axes2, ' all ' );plot2 = plot(nu,abs(eval(funcStr) - data2plot), ' Parent ' ,axes2);set(plot2(1), ' LineWidth ' ,1, ' Color ' ,[0 0 0]);xlabel( ' Parameter \it{\nu} ' , ' FontSize ' ,14);ylabel( ' Absolute difference ' , ' FontSize ' ,14);endend Conclusion
In this work we derived a rational approximation of the Fourier transformthat with help of a CAS can be readily rearranged as
F { f ( t ) } ( ν ) ≈ P ( ν ) Q ( ν ) . This method of the rational approximation is based on integration involv-ing exponential decay multiplier e − σt . The computational test we performedshows that this method of the Fourier transform can provide relatively ac-curate approximations of the Fourier transform even for the functions withdiscontinuities like rect ( t ) and it rect ( t ). Furthermore, this method showsthat for the well-behaved function f ( t ) = √ πe − ( πt ) + i (cid:104) π / te − ( πt ) (cid:105) withonly 16 summation terms the rational approximations (14) and (15) providethe Fourier transform with absolute differences not exceeding 3 × − and9 × − for its real and imaginary parts, respectively. Acknowledgments
This work is supported by National Research Council Canada, Thoth Tech-nology Inc., York University and Epic College of Technology.
References [1] M. Abramowitz and I.A. Stegun. Error function and Fresnel integrals.Handbook of mathematical functions with formulas, graphs, and mathe-matical tables. 9th ed. New York 1972, 297-309.[2] S.M. Abrarov and B.M. Quine, Sampling by incomplete cosine expan-sion of the sinc function: Application to the Voigt/complex error func-tion, Appl. Math. Comput., 258 (2015) 425-435. http://dx.doi.org/10.1016/j.amc.2015.01.072 [3] S.M. Abrarov and B.M. Quine, A rational approximation for efficientcomputation of the Voigt function in quantitative spectroscopy, J. Math.Res., 7(2) (2015) 163-174. http://dx.doi.org/10.5539/jmr.v7n2 https://doi.org/10.1016/j.apnum.2019.08.030 [5] B.H. Armstrong, Spectrum line profiles: The Voigt function, J. Quant.Spectrosc. Radiat. Transfer, 7(1) (1967) 61-88. https://doi.org/10.1016/0022-4073(67)90057-X [6] A. Berk, Voigt equivalent widths and spectral-bin single-line transmit-tances: Exact expansions and the MODTRAN ® https://doi.org/10.1016/j.jqsrt.2012.11.026 [7] R.N. Bracewell, The Fourier transform and its applications, 3rd ed.,McGraw-Hill, 2000.[8] G. Colldeforns-Papiol and L. Ortiz-Gracia, Computation of market riskmeasures with stochastic liquidity horizon, J. Comput. Appl. Math., 342C(2018) 431-450, https://doi.org/10.1016/j.cam.2018.03.038 [9] W.B. Gearhart and H.S. Shultz, The function sin ( x ) /x , College Math. J,21 (1990) 90-99. [10] E.W. Hansen, Fourier transforms. Principles and applications, John Wi-ley & Sons, 2014.[11] M. Kac, Statistical independence in probability, analysis and numbertheory. Washington, DC: Math. Assoc. Amer., 1959.[12] S. C. Maree, L. Ortiz-Gracia and C. W. Oosterlee, Fourier and waveletoption pricing methods. High-performance computing in finance, prob-lems, methods, and solutions, edited by M. A. H. Dempster et al. , CHAP-MAN & HALL/CRC 2018, pp. 249-272.[13] L. Ortiz-Gracia and C. W. Oosterlee, SIAM J. Sci. Comput., A highlyefficient Shannon wavelet inverse Fourier technique for pricing Euro-pean options, 38 (1) (2016) B118B143. https://doi.org/10.1137/15M1014164 [14] G.B. Rybicki, Dawson’s integral and the sampling theorem, Comp.Phys., 3 (1989) 85-87. http://dx.doi.org/10.1063/1.4822832 http://dx.doi.org/10.1016/j.jqsrt.2013.04.020http://dx.doi.org/10.1016/j.jqsrt.2013.04.020