Peaks and jumps reconstruction with B-splines scaling functions

# Peaks and jumps reconstruction with B-splines scaling functions

Luis Ortiz-Gracia Centre de Recerca Matemàtica, Campus de Bellaterra, Edifici C, 08193 Bellaterra (Barcelona), Spain  and  Josep J. Masdemont Departament de Matemàtica Aplicada I, Universitat Politècnica de Catalunya, Diagonal 647, 08028 Barcelona, Spain
October 2012
###### Abstract.

We consider a methodology based in B-splines scaling functions to numerically invert Fourier or Laplace transforms of functions in the space . The original function is approximated by a finite combination of order B-splines basis functions and we provide analytical expressions for the recovered coefficients. The methodology is particularly well suited when the original function or its derivatives present peaks or jumps due to discontinuities in the domain. We will show in the numerical experiments the robustness and accuracy of the method.

## 1. Introduction

Fourier and Laplace transforms are useful in a wide number of applications in science and engineering. As it is well known, the Fourier transform is closely related to the Laplace transform for zero value funcions on the negative time axis. There is a strong interest in the efficient numerical inversion of Laplace transforms ([1],[2]) and Fourier transforms, due to the fact that the solutions to some problems are known in the transform domain rather than in the original domain. An increasing number of papers have recently appeared to invert Fourier transforms with wavelets, like [8] and [9] with coiflets wavelets and [10] with Mexican, Morlet, Poisson and Battle-Lemarié wavelets. In particular in the Financial Engineering context, [11] inverts a Laplace transform by means of B-splines wavelets of order .

Recently, a new method called the COS method developed in [7] for solving the inverse Fourier integral is capable to accurately recover a function from its Fourier transform in a short CPU time being as well of very easy implementation. However, when the function to be recovered presents discontinuities or it is highly peaked, a lot of terms in the expansion must be considered to reduce the approximation error.

In this paper we present a novel approximation based on B-splines scaling functions to numerically invert Fourier transforms that it is particularly suited for functions that exhibit peaks and/or jumps in its domain. There are some properties that, at first glance, make these basis functions particularly suited to approximate such non-smooth functions. B-splines are the most regular scaling functions with the shortest support for a given polynomial degree. Another important fact is the explicit formulation in the time (or space) domain as well as in the frequency domain. However a slight drawback of these basis functions is that the system that they form is semi-orthogonal, i.e., the scaling functions are orthogonal among different scales but not necessarily at the same scale.

In previous work [14] the authors numerically inverted the Laplace transform of a distribution function in the interval making use of the Haar scaling functions. Now, in the present paper, we consider the problem of inverting the Fourier transform to recover an function by approximating it by a finite sum of B-splines scaling functions of order (where is the particular case of the Haar system). We also provide a list of the different errors accumulating within the numerical procedure. So, here, the Fourier inversion is carried out with B-splines scaling functions rather than with wavelet functions ([10]) or a combination of both ([8],[9]). We fix the scale parameter in the wavelet expansion and only remaining is the translation parameter, facilitating the inversion procedure. Furthermore, we provide an analytical expression for the coefficients of the approximation.

As will be shown in the section devoted to numerical examples, the Wavelet Approximation (WA) that we present is well capable to detect jumps or peaks produced by discontinuities in the function itself or in first derivatives. On the contrary, the COS method is better to approximate analytical functions.

The paper is organized as follows. In Section 1.1 we give a short literature review regarding the Laplace transform inversion. Section 2 gives a brief introduction concerning multiresolution analysis and B-splines scaling functions. In Section 3 we present the Wavelet Approximation method to recover functions from its Fourier transform by means of B-splines scaling functions. Section 4 is devoted to the COS method, while numerical experiments, comparing the Wavelet Approximation and the COS method are shown in Section 5. Finally, Section 6 concludes.

### 1.1. Laplace transform inversion

Suppose that is a real- or complex-valued function of the variable and is a real or complex parameter. We define the Laplace transform of as,

 (1) ˜f(s)=∫+∞0e−sxf(x)dx=limτ→+∞∫τ0e−sxf(x)dx,

whenever the limit exists (as a finite number).

We state the inverse transform as a theorem (see [6] for a detailed proof).

###### Theorem 1.

(Bromwich inversion integral) If the Laplace transform of exists, then,

 (2) f(x)=limk→+∞12πi∫σ+ikσ−ik˜f(s)esxds,x>0,

where for some positive real number and is any other real number such that .

The usual way of evaluating this integral is via the residues method taking a closed contour, often called the Bromwich contour.

In this section we present some numerical algorithms to invert the Laplace transform. A natural starting point for the numerical inversion of Laplace transforms is the Bromwich inversion integral stated in Theorem 1. If we choose a specific contour and perform the change of variables in (2), we obtain an integral of a real valued function of a real variable. Then, after algebraic manipulation and after applying the Trapezoidal Rule we obtain,

 (3) f(x)≃fh(x)=heσxπR(˜f(σ))+2heσxπ+∞∑k=1R(˜f(σ+ikh))cos(khx),

where denotes the real part of . A detailed analysis of the errors can be found in [2].

As pointed out in [2], the Bromwich inversion integral is not the only inversion formula and there are quite different numerical inversion algorithms. We refer the readers to the Laguerre series representation given in [1], which is known to be an efficient method for smooth functions. However, if is not smooth at just one value of , then the Laguerre method has difficulties at any value of .

In the context of numerical Laplace inversion, [8] recovers the function with a procedure based on wavelets. They consider in expression (2), where is a real variable and is a real constant that fulfills , assuming that when . Then, equation (1) can be rewritten as,

 ˜f(β+iω)=∫+∞−∞e−βxe−iωxf(x)dx.

Defining,

 (4) h(x)=f(x)e−βx,then^h(ω)=˜f(β+iω),

where denotes the Fourier transform of . The authors expand in terms of Coiflets wavelets,

 (5) ^h(ω)=+∞∑k=−∞cm,kϕm,k(ω)++∞∑j=m+∞∑k=−∞dj,kψj,k(ω).

where , , being and the scaling and wavelet functions respectively.

The next step consists in inverting the expression (5) by means of the Fourier inversion formula. Finally, considering the expression (4), the formulae of Laplace inversion become,

 fm(x)=eβx2m+1π^ϕ(−x2m)+∞∑k=−∞˜f(β+iM1+k2m)eixk/2m,f(x)=limm→+∞fm(x).

One drawback of this approximation is that the wavelet approach involves an infinite product of complex series and the computation of the Fourier transform of some scaling functions. This can look intimidating for practical applications and may also take relatively long computational time.

Based on operational matrices and Haar wavelets, the author in [4] presents a new method for performing numerical inversion of the Laplace transform where only matrix multiplications and ordinary algebraic operations are involved. However, the essential step in the method consists in expressing the Laplace transform in terms of , which is impossible when we just know numerically the transform. Another drawback of this method is that the matrices become very large for larger scales.

## 2. Multiresolution analysis and cardinal B-splines

A natural and convenient way to introduce wavelets is following the notion of multiresolution analysis (MRA). Here we provide the basic definitions and properties regarding MRA and B-spline wavelets, for further information see [13, 5].

###### Definition 1.

A countable set of a Hilbert space is a Riesz basis if every element of the space can be uniquely written as , and there exist positive constants and such that,

 A∥f∥2≤∑n|cn|2≤B∥f∥2.
###### Definition 2.

A function is called a scaling function, if the subspaces of , defined by,

 Vm=closL2(R){ϕm,k:k∈Z},m∈Z,

where , satisfy the properties,

1. .

2. .

3. .

4. For each , is a Riesz (or unconditional) basis of .

We also say that the scaling function generates a multiresolution analysis of .

The order cardinal -spline function, , is defined recursively by a convolution:

 Nj(x)=∫∞−∞Nj−1(x−t)N0(t)dt=∫10Nj−1(x−t)dt,j≥1,

where,

 N0(x)=χ[0,1)(x)={1if x∈[0,1)0otherwise.

Alternatively,

 Nj(x)=xjNj−1(x)+j+1−xjNj−1(x−1),j≥1.

We note that cardinal -spline functions are compactly supported, since the support of the order -spline function is , and they have as the Fourier transform,

 ˆNj(w)=(1−e−iwiw)j+1.

In this paper we consider as the scaling function which generates a MRA (see Figure 1). Clearly, for we have the scaling function of the Haar wavelet system. We also remark that from the previous discussions, for every function , there exists a unique sequence , such that,

 fm(x)=∑k∈Zcjm,kϕjm,k(2mx−k).

## 3. The wavelet approximation method

Let us now consider a function and its Fourier transform, whenever it exists:

 ˆf(w)=∫+∞−∞e−iwxf(x)dx.

Since we can expect that decays to zero, so it can be well approximated in a finite interval by,

 (6) fc(x)={f(x)if x∈[a,b],0otherwise.

Let us approximate for all , where,

 (7) fcm,j(x)=(j+1)⋅(2m−1)∑k=0cjm,kϕjm,k((j+1)⋅x−ab−a),j≥0,

with convergence in -norm. Note that we are not considering the left and right boundary scaling functions (we refer the reader to Section 3 in [12] for a detailed description of scaling functions on a bounded interval).

The main idea behind the Wavelet Approximation method is to approximate by and then to compute the coefficients by inverting the Fourier Transform. Proceeding this way,

Introducing the change of variables , gives us,

 ˆf(w)≃b−aj+1⋅e−iaw(j+1)⋅(2m−1)∑k=0cjm,k∫+∞−∞e−iwb−aj+1yϕjm,k(y)dy=b−aj+1⋅e−iaw(j+1)⋅(2m−1)∑k=0cjm,kˆϕjm,k(b−aj+1⋅w).

Finally, taking into account that and performing the change of variables , we have,

 ˆf(2m(j+1)b−ai⋅log(z))≃2−m2b−aj+1⋅z2m(j+1)ab−aˆϕj(i⋅log(z))(j+1)⋅(2m−1)∑k=0cjm,kzk.

If we consider,

 Pm,j(z)=(j+1)⋅(2m−1)∑k=0cjm,kzkandQm,j(z)=2m2(j+1)z−2m(j+1)ab−aˆf(2m(j+1)b−ai⋅log(z))(b−a)ˆϕj(i⋅log(z)),

then, according to the previous formula, we have,

 (8) Pm,j(z)≃Qm,j(z).

Since is a polynomial, it is (in particular) analytic inside a disc of the complex plane for . We can obtain expressions for the coefficients by means of the Cauchy’s integral formula. This is,

 cjm,k =12πi∫γPm,j(z)zk+1dz,k=0,...,(j+1)⋅(2m−1),

Considering now the change of variables , , gives us,

 (9) cjm,k=12πrk∫2π0Pm,j(reiu)eikudu=12πrk∫2π0[R(Pm,j(reiu))cos(ku)+I(Pm,j(reiu))sin(ku)]du,

where , and and stand for the real and imaginary parts of , respectively.

Note that if then we can further expand the expression above by considering,

 (10) cjm,k=2πrk∫π0R(Pm,j(reiu))cos(ku)du.

On the other side, since , we have,

 (11) Qm,j(z)=2m2(j+1)z−2m(j+1)ab−aˆf(2m(j+1)b−ai⋅log(z))(log(z))j+1(b−a)(z−1)j+1,

and it has a pole at . Finally, making use of (8) and taking into account the former observation, we can exchange by in (9) and (10) to obtain, respectively,

 (12) cjm,0≃12π∫2π0R(Qm,j(reiu))du,

and,

 (13) cjm,k ≃2πrk∫π0R(Qm,j(reiu))cos(ku)du,k=1,...,(j+1)⋅(2m−1),

where is a positive real number.

In practice, both integrals in (12) and (13) are computed by means of the Trapezoidal Rule, and we can define,

 I(k) =∫π0R(Qm,j(reiu))cos(ku)du, I(k;h) =h2(Qm,j(r)+(−1)kQm,j(−r)+2M−1∑j=1R(Qm,j(reihs))cos(khs)),

where and for all . Proceeding this way we find,

 (14) cjm,k≃2πrkI(k)≃2πrkI(k;h)=1Mrk(Qm,j(r)+(−1)kQm,j(−r)+2M−1∑s=1R(Qm,j(reihs))cos(khs)),

where .

Let us summarize four sources of error in our procedure to compute the numerical Fourier transform inversion using cardinal B-splines wavelets. These are:

1. Truncation of the integration range,

 E1(x):=f(x)−fc(x),x∈R∖[a,b].
2. The approximation error at scale ,

 E2(x):=fc(x)−fcm,j(x),x∈[a,b].
3. The discretization error, which results when approximating the integral by using the trapezoidal rule. We can apply the formula for the error of the compound trapezoidal rule considering,

 qjm,k(u)=R(Qm,j(reiu))cos(ku),E3:=I(k)−I(k;h),

and assuming that . Then,

 (15) |E3|=π312M2∣∣(qjm,k(μ))′′∣∣,μ∈(0,π).
4. The roundoff error. If we can calculate the sum in expression (14) with a precision of , then the roundoff error after multiplying by a factor is approximately . Then, the roundoff error increases when approaches to 0.

### 3.1. Choice of the parameter r

As mentioned before, the choice of the parameter may have influence in both, the discretization and the roundoff errors. In this section we present a detailed analysis of the errors listed before in order to determine the optimum value. To do this, let us consider a step function defined in , where represents the indicator function, and its Fourier transform .

Due to the shape of it seems that the best B-splines basis to perform the approximation is based in the Haar scaling functions (B-splines of order ). In this particular case and its derivatives up to order can be computed relatively straightforwardly, so that we will be able to calculate the optimal value of the parameter in order to minimize the discretization and the roundoff errors. We also demonstrate that the approximation error (type (B)) is in this case. To do so, we consider,

 Qm,0(reiu)=r2m(cos(2mu)+isin(2mu))−r2m−1(cos(2m−1u)+isin(2m−1u))2m/2(rcos(u)+irsin(u)−1),

where , and,

 R(Qm,0(reiu))==r2m+1cos(2mu−u)−r2mcos(2mu)−r2m−1+1cos(2m−1u−u)+r2m−1cos(2m−1u)2m/2(r2−2rcos(u)+1).

Now, we must choose an appropriate value to control both the discretization and the round-off errors. First of all, we consider the discretization error which can be estimated by means of expression (15). We note that since,

 0<(r−1)2≤r2−2rcos(u)+1≤(r+1)2,∀u∈[0,π],r≠1.

So we have,

 dduq0m,k(u)=dduR(Qm,0(reiu))cos(ku)−kR(Qm,0(reiu))sin(ku),d2du2q0m,k(u)=d2du2R(Qm,0(reiu))cos(ku)−2kdduR(Qm,0(reiu))sin(ku)−k2R(Qm,0(reiu))cos(ku).

and,

 (16) ∣∣R(Qm,0(reiu)∣∣≤r2m+1+r2m+r2m−1+1+r2m−12m/2(r−1)2≃r2m−1+o(r2m−1),

where the last approximation holds for suitably small values of the parameter . For sake of simplicity, we consider only the terms with smaller exponents in the parameter for the expressions and . Then,

 (17) ∣∣∣dduR(Qm,0(reiu))∣∣∣≤2(m−2)/2r2m−1+A(r)(r−1)4≃r2m−1+o(r2m−1),

and,

 (18) ∣∣∣d2du2R(Qm,0(reiu))∣∣∣≤2(3m−4)/2r2m−1+B(r)(r−1)8≃r2m−1+o(r2m−1),

where and are polynomials in with degree greater than , and the approximations in (17) and (18) hold for suitably small values of the parameter . Finally, taking into account expressions (15),(16),(17) and (18) we have,

 |E3|≤π312M2(r2m−1+2kr2m−1+k2r2m−1)+o(r2m−1)=π312M2(k2+2k+1)r2m−1+o(r2m−1).

We note that as while the roundoff error is increasing in as approaches to zero.

Then, the total error should be approximately minimized when the two estimates are equal. This leads to the equation,

After algebraic manipulation, we find,

 (19) rm,k=(12M⋅10−ηπ3(k2+2k+1))12m−1+k,k=0,…,2m−1.

As mentioned before, the roundoff error arises when multiplying the sum in expression (14) by the pre-factor . Let us consider , then the of interest is which is the greatest value that this parameter can take (small values of do not cause roundoff errors).

The left plot of Figure 2 represents the pre-factor for values of , while the right plot shows the pre-factor values for . We also display in Table 1 the pre-factor values for different values of and scales , and .

On the one hand, it is clear that we must concentrate in this second interval in order to get a reasonable roundoff error and on the other hand, the discretization error grows when is close to . Later, in the numerical examples section, we will confirm the theory developed above for the step function.

Since is compactly supported, we have .

Note that in the case that ,

where, is the Haar basis (orthonormal) system in . Then,

 (20) ∥E2∥2L2([a,b])=∥fc−fcm∥2L2([a,b])=∥∥ ∥∥+∞∑l=m2l−1∑k=0d0l,kψl,k(x−ab−a)∥∥ ∥∥2L2([a,b])=(b−a)⋅+∞∑l=m2l−1∑k=0|d0l,k|2,

since . Then, the approximation error depends on the length of the interval and the detail coefficients,

 (21) d0l,k=∫Rfc(x)⋅ψl,k(x)dx.

Furthermore, since the detail coefficients (21) are zero, then we have also that and the approximation is exact at any scale level .

## 4. The COS method

For completeness, we present here the methodology developed in [7] for solving the inverse Fourier integral111In order to maintain the notation used by the authors, here (22) represents the characteristic function, and hence the Fourier transform of a density function .. The main idea is to reconstruct the whole integral from its Fourier-cosine series expansion extracting the series coefficients directly from the integrand. Fourier-cosine series expansions usually give an optimal approximation of functions with a finite support [3]. In fact, the cosine expansion of in equals the Chebyshev series expansion of in .

For a function supported on , the cosine expansion reads,

 f(θ)=A02++∞∑k=1Akcos(kθ),

with . For functions supported on any other finite interval , the Fourier-cosine series expansion can easily be obtained via a change of variables,

 θ:=x−ab−aπ,x=b−aπθ+a.

 f(x)=A02++∞∑k=1Akcos(kπx−ab−a),

with,

 (23) Ak=2b−a∫baf(x)cos(kπx−ab−a)dx.

Since any real function has a cosine expansion when it is finitely supported, the derivation starts with a truncation of the infinite integration range in (22). Due to the conditions for the existence of a Fourier transform, the integrands in (22) have to decay to zero at and we can truncate the integration range in a proper way without losing accuracy.

Suppose is chosen such that the truncated integral approximates the infinite counterpart very well, i.e.,

 (24) ξ1(w):=∫baeiwxf(x)dx≃∫Reiwxf(x)dx=ξ(w).

Here, denotes a numerical approximation.

Comparing equation (24) with the cosine series coefficients of on in (23), we find that,

 Ak≡2b−aR(ϕ1(kπb−a)e−ikaπb−a),

where denotes the real part of the argument. It then follows from (24) that with,

 Fk≡2b−aR(ϕ(kπb−a)e−ikaπb−a).

We now replace by in the series expansion of on , i.e.,

 (25) f1(x)=A02++∞∑k=1Fkcos(kπx−ab−a),

and truncate the series summation such that,

 (26) f2(x)=A02+N−1∑k=1Fkcos(kπx−ab−a).

The resulting error in consists of two parts, a series truncation error from (25) to (26) and an error originated from the approximation of by . An error analysis that takes these different approximations into account is presented in [7].

The COS method performs well when approximating smooth functions, but many terms are needed in case that the function or its first derivative has discontinuities along the domain of approximation.

## 5. Numerical examples

The aim here is to show the accuracy of B-splines to invert Fourier transforms of extremely peaked or discontinuous functions with finite support. We will consider the B-splines scaling functions of order . We denote by WAi-j the Wavelet Approximation method with B-splines of order at scale , which involves coefficients and COS-N the COS approximation method with terms.

To compute the coefficients of the Fourier inversion in expression (14), we must set the parameters. For this purpose we consider , where is the order of the B-spline considered and is the scale parameter. Observe that if we take instead of in (14) this leads to the following expression,

 (27) cm,k≃12krk(Qm,j(r)+(−1)kQm,j(−r)+22k−1∑s=1R(Qm,j(reisπ2k))cos(ksπ2k))=12krk(Qm,j(r)+(−1)kQm,j(−r)+2k−1∑s=1(−1)sR(Qm,j(reisπ2k))),

for , so the computation time reduces for large scale approximations.

Finally, we set . Although we know that and , we must take into account the two types of errors (C) (discretization) and (D) (roundoff) listed previously, which may have influence on the computation of the coefficients. Due to the fact that the function is in general intractable from an analytical point of view, we carried out intensive simulations in order to understand the influence of parameter . We did these simulations considering the test functions and at different scale levels and used B-splines scaling functions of order . To show just an example, we have plotted the results for and B-splines of order at scale (Figure 3) and (Figure 4). The colors represent the magnitude of the logarithm of the absolute errors. As we can observe, the absolute error remains constant for values . When , the error increases for high values of (i.e. high values of the translation parameter ) and decreases for low values of (i.e. low values of the translation parameter ). It is worth mentioning that for high scales, the error grows very rapidly when (empirical findings demonstrate that , for this reason we take approximately the midpoint of the interval avoiding ), while this effect diminishes for shorter scales. In fact, as we will see later with the step function, the roundoff error almost disappears at very low scales (for instance ) and the discretization error tends to zero when tends to zero. These facts confirm the high impact of the roundoff error at high scale levels and only the impact of the discretization error at very low scales.

For comparison, we gave in Section 4 a brief overview of the COS method developed by Fang and Oosterlee ([7]) that is based on a Fourier-cosine series expansion and which usually gives optimal approximations of functions with finite support ([3]). The COS method is state of the art and has been applied to efficiently recover density functions from their Fourier transforms in order to solve important problems arising in Computational Finance.

Let us consider the step function presented before. Observe that if we take , there is almost no roundoff error and the only remaining error is the discretization error . If we assume that , then according to (19) , . The left plot of Figure 5 shows the approximation to the step function with the COS and the WA method with . The right plot of Figure 5 represents the absolute error of the approximation with the COS and the WA method with and the optimum computed with expression (19). Note that with just coefficients the WA method is capable to accurately approximate the step function, while the COS method suffers from the Gibbs phenomenon even if a lot of terms (we show up to ) are added to the COS expansion.

Note that if we consider the approximation to the step function at scale , then the values computed for the parameter are all of them in a neighborhood of in accordance with the massive simulations performed before.

The conclusion in that case is that the COS method performs poorly around the jump and exhibits the Gibbs phenomenon, while Haar wavelets are naturally capable to deal with this discontinuity.

### 5.1. Exponential function

Let us consider the exponential function, and its Fourier transform .

Observe that this function becomes extremely peaked when . We will test the inversion with and .

We consider the interval and . We carry out the comparison between both methods using approximately the same number of terms in the approximation. It is worth mentioning that the computation of wavelet coefficients is more time consuming than the calculation of COS coefficients. Moreover, the COS method is easier to implement. Results are reported in Table 2 and plotted in Figure 6. We observe that linear B-splines are the most suitable basis functions to approximate highly peaked functions like the exponential we are considering. While adding many more terms in the COS expansion improves only a bit the approximation, when we consider B-splines of order at higher scales the approximation error decays much more faster. The WA method with B-splines of order performs better than the WA method with B-splines of orders and .