### High Accuracy Reconstruction from Wavelet Coefficients

Fritz Keinert* and Soon-Geol Kwon

*Department of Mathematics, Iowa State University, Ames, Iowa 50011*

*Communicated by Gregory Beylkin*

Received June 10, 1996; revised February 28, 1997

*The accuracy of the wavelet approximation at resolution h Å 2** ^{0n}*to a smooth

*function f is limited by O(h*

^{N}*), where N is the number of vanishing moments of*

*the mother wavelet c. For any positive integer p, we derive a new approximation*

*formula which allows us to recover a smooth f from its wavelet coefficients*

*with accuracy O(h*

^{p}*). Related formulas for recovering derivatives of f are also*given. q1997 Academic Press

**1. INTRODUCTION**

Let f, f_{H} be the scaling functions of a biorthogonal multiresolution approximation.

The projection

*P*H^{n}*f* Å

### ∑

^{`}

*jÅ0*`

»*f, f**n,j*…f_{H}*n,j* (1.1)

*of a function f*√*L*^{2}(R*) onto the space V˜** _{n}* spanned by

f_{H}*n,k**(x)* Å2* ^{n/2}*f

_{H}(2

^{n}*x*0

*k),*

*k*√Z, (1.2)

*is interpreted as an approximation to f at resolution h*Å 2

*. Such projections form the starting point and/or result of any numerical method based on wavelets, from signal processing to wavelet – Galerkin methods.*

^{0n}*In practice, the connection between point values of the smooth original function f*
and its scaling function coefficients is often established by simple formulas such as

*f**n,j*Å »*f, f**n,j*… É2^{0n/2}*f ((j*/*t)h),*

* E-mail: keinert@iastate.edu.

293

where t Å

### *

*xf(x)dx is the center of mass of f. For many applications, this is*perfectly adequate.

*Higher accuracy formulas for the computation of f**n,j*have been published (see, e.g.,
*[6]), but the problem of recovering accurate point values of f from its coefficients f**n,j*

has not received much attention.

*It is well known that if the wavelet c has N vanishing moments and f belongs to*
*C** ^{N}*, then

É*f(x)*0*P*H^{n}*f (x)*^{É} Å*O(h** ^{N}*).

In some applications, such as wavelet – Galerkin methods, we may know the coeffi-
cients of a smooth solution function to high accuracy. How do we recover point values
*of f, and possibly some derivatives of f, to comparable accuracy?*

*For any given scaling function f and accuracy p, we find in this paper a function*
*b with support in an interval of length p so that for smooth f,*

É*f (x)*0

### ∑

*j*

»*f, f**n,j*…b*n,j**(x)*^{É} Å*O(h** ^{p}*), (1.3)

where b*n,j*Å2* ^{n/2}*b(2

^{n}*x*0

*k), k*√Z, in analogy to (1.2). Related functions b

*can be*

^{[r]}*constructed to recover the rth derivative of f to accuracy O(h*

^{p0r}*) for any r*õ

*p.*

While the series in (1.3) looks similar to a scaling function expansion, b and b^{[r]}

do not satisfy any recursion relations in general.

Numerical experiments indicate that our alternative reconstruction formulas perform
*significantly better than the standard scaling function series (1.1) for smooth f and*
*are no worse for highly oscillating or nonsmooth f. In addition, the new approach*
*allows the reconstruction of derivatives of f.*

**2. WAVELET BASICS**

The fundamentals of wavelet theory have been described in many publications, so we refer to the literature for most of the background, and we only introduce the specific notation and results we need. A standard reference for wavelets is [3]. Biortho- gonal wavelets are discussed in more detail in [2]. Throughout this paper, we will be using biorthogonal wavelets, which include orthogonal wavelets as a special case.

*L*^{2}(_{R}) is the space of square integrable complex-valued functions on_{R}with inner
product

»*f, g*… Å

### *

^{R}

*f (x)g(x)dx,*

*where the bar denotes complex conjugation. l*2 is the space of square summable
complex-valued sequences with inner product

»*x, y*… Å

### ∑

^{`}

*iÅ0*`

*x**i**y**i*.

*C*^{r}*is the space of r times continuously differentiable functions on* _{R}.

*A multiresolution analysis (MRA) [5] of L*^{2}(R*) is a sequence {V**j*}*j√*Z of closed
*subspaces of L*^{2}(R) with the properties

• *V**j*, *V**j/1*

• <^{`}*jÅ0`* *V**j**is dense in L*^{2}(R), and>^{`}*jÅ0`* *V**j*Å{0}

• *f (x)*√*V**j*B*f (2x)*√*V**j/1*

• *f (x)*√*V**j*B*f (x*02^{0j}*k)* √*V**j**for all k*√Z

• *there exists a scaling function f(x)* √*V*0 *such that the set of translates {f(x*0
*k)}**k√*Z *forms a Riesz basis of V*0.

*Let W*0 *be a complement of V*0*in V*1, that is,
*V*1Å*V*0/*W*0.

*A function c(x)* √ *W*0 *whose translates {c(x*0 *k)}**k√*Z *form a Riesz basis of W*0is
*called a mother wavelet. We assume that such a function exists.*

*Since V*0*, W*0*are both subspaces of V*1*, recursion relations*
*f(x)* Å

q

2

### ∑

*k*

*h**k**f(2x* 0*k),*

*c(x)* Å^{q}2

### ∑

*k*

*g**k**f(2x* 0*k),* (2.1)

*must hold for some sequences {h**k**}, {g**k**} in l*2.
We denote the scaled translates of f, c by

f*n,k**(x)*Å2* ^{n/2}*f(2

^{n}*x*0

*k),*

c*n,k**(x)*Å2* ^{n/2}*c(2

^{n}*x*0

*k).*(2.2)

*Biorthogonal wavelets are defined by two MRAs of L*

^{2}(R

*), denoted by {V*

*j*} and

*{V˜*

_{j}*}, whose basis functions satisfy the biorthogonality relations*

»f*n,j*, f_{H} *n,k*… Åd*jk*,

»c*n,j*, c_{H} *n,k*… Åd*jk*,

»f*n,j*, c_{H}*n,k*… Å »c*n,j*, f_{H}*n,k*… Å 0. (2.3)
Here d*jk* is the Kronecker delta

d*jk*Å

### H

^{0}

^{1}

^{if j}^{if j}^{x}

^{Å}

^{k,}^{k.}*For any function f*√*L*^{2}(R*), projections P**n**f, Q**n**f, P˜*_{n}*f, Q˜*_{n}*f of f onto V**n**, W**n**, V˜*_{n}*, W˜** _{n}*,
respectively, are given by

*P**n**f*Å

### ∑

^{`}

*jÅ0*`

»*f, f*_{H}*n,j*…f*n,j*, *P*H^{n}*f*Å

### ∑

^{`}

*jÅ0*`

»*f, f**n,j*…f_{H}*n,j*,

*Q**n**f*Å

### ∑

^{`}

*jÅ0*`

»*f, c*_{H}*n,j*…c*n,j*, *Q*H^{n}*f* Å

### ∑

^{`}

*jÅ0*`

»*f, c**n,j*…c_{H}*n,j*.

*In many applications, P**n**f and P˜*_{n}*f are interpreted as approximations to f at resolution*
*h*Å2^{0n}*, while Q**n**f and Q˜*_{n}*f represent the fine detail in f at resolution h.*

*The relationship between {V**n**}, {V˜** _{n}*} etc. is completely symmetric. To any definition
or algorithm there is a dual, with the tildes in opposite places. We choose to work

*with P˜*

_{n}*instead of P*

*n*in the following sections because it keeps the notation simpler.

*The moments of f, c, f*_{H}, c_{H} are defined as

*M**i*Å

### *

^{x}

^{i}

^{f(x)dx,}

^{M}^{H}

^{i}^{Å}

### *

^{x}

^{i}^{f}

^{H}

^{(x)dx,}*N**i*Å

### *

^{x}

^{i}

^{c(x)dx,}

^{N}^{H}

^{i}^{Å}

### *

^{x}

^{i}^{c}

^{H}

^{(x)dx.}^{(2.4)}

We assume that the scaling functions have been normalized so that*M*0Å *M*H ^{0}^{Å} 1.

*We say that c has N vanishing moments ifN*0Å*N*1Å rrr Å*N**N01*Å0,*N**N* x0.

*For later use, we also define the shifted moments of f,*

*M**i,j*Å

### *

^{x}

^{i}

^{f(x}^{0}

^{j)dx}^{Å}

### *

^{(x}^{/}

^{j)}

^{i}

^{f(x)dx}^{Å}

*sÅ0*

^{∑}

^{i}### S

^{i}^{s}### D

^{j}

^{s}

^{M}

^{i0s}^{.}

^{(2.5)}

It is easy to verify that

*M**k,j/1*Å

### ∑

^{k}*lÅ0*

### S

^{k}^{l}### D

^{M}

^{l,j}^{.}

^{(2.6)}

**3. CONSTRUCTION OF ALTERNATIVE RECONSTRUCTION FORMULAS**
*Assume that for some function f and for some biorthogonal wavelet basis we know*
*the projection P˜*_{n}*f onto V˜** _{n}*,

*P*H^{n}*f*Å

### ∑

^{`}

*jÅ0*`

*f**n,j*f_{H}*n,j*, (3.1)

where

*f**n,j*Å »*f, f**n,j*….

*How well can we recover the point values of f from this?*

*Assume that c has N vanishing moments. If f* √*C** ^{N}*, it is well known and easy to

*prove that for any point x,*

É*f (x)*0 *P*H^{n}*f (x)*^{É} Å*O(h** ^{N}*),

*where h*Å2^{0n}*. (See, e.g., [6, Eq. (1.11)]; note that the meaning of N, N˜ is reversed,*
compared to our paper).

*There are more detailed results available about the convergence rate of P˜*_{n}*f to f*
*under various conditions of f and f (see, for example, Walter [7], Kelly et al. [4]),*
but they are all based on the original scaling function series (3.1).

We propose instead to use a different series,

*B**n**f*Å

### ∑

^{`}

*jÅ0*`

*f**n,j*b*n,j*, (3.2)

where

b*n,j**(x)* Å2* ^{n/2}*b(2

^{n}*x*0

*j)*(3.3)

*in analogy to (2.2). For any choice of f and for any p*√N, we will produce a function

*b(x) with support of length p so that for f*√

*C*

*,*

^{p}É*f (x)*0*B**n**f (x)*^{É} Å*O(h** ^{p}*). (3.4)

*For simplicity, the dependence of b on p and f is not usually expressed in the notation.*

*If necessary we will write b(x; p) and, similarly, for dependence on other parameters.*

More generally, we want to find functions b* ^{[r]}*so that

É*f*^{(r)}*(x)* 0*B*^{[r]}*n* *f (x)*^{É} Å*O(h** ^{p0r}*), (3.5)
where

*f** ^{(r)}* Å

*d*

^{r}*dx*

^{r}*f,*

*B*

^{[r]}*n*

*f*Å

### ∑

^{`}

*jÅ0*`

*f**n,j*b^{[r]}*n,j*,

b^{[r]}*n,j**(x)* Å2* ^{n/2}*b

*(2*

^{[r]}

^{n}*x*0

*j).*(3.6)

*To construct b, fix a scaling function f and integers n, p and assume that we are*

*given the scaling function coefficients f*

*n,j*

*, j*Å

*0, . . . , p*01. Following a standard

*approach from numerical analysis, we first attempt to find coefficients c*

*n,j*

*(x), j*Å0,

*. . . , p*01, so that

*f (x)*Å^{p01}

### ∑

*jÅ0*

*f**n,j**c**n,j**(x)* *for f (x)*Å*x*^{k}*, k*Å*0, . . . , p*01. (3.7)

*For f (x)*Å*x** ^{k}*,

*f**n,j*Å »*f, f**n,j*… Å

### *

^{x}

^{k}^{2}

^{n/2}^{f(2}

^{n}

^{x}^{0}

^{j)dx}^{Å}

^{h}

^{k/(1/2)}

^{M}

^{k,j}^{.}

Equation (3.7) leads to a system of linear equations,

**H***n***Mc***n**(x)* Å* p(x),* (3.8)

where

**H***n*Å **H***n**(p)* Å*h*^{1/2}
*h*^{0}

0
*h*^{1}

???

0

*h** ^{p01}*
,

**M**Å* M(f, p)* Å

*M*0,0 *M*0,1 ??? *M**0,p01*

*M*1,0 *M*1,1 ??? *M**1,p01*

: : ??? :

*M**p01,0* *M**p01,1* ??? *M**p01,p01*

,

**c***n**(x)*Å**c***n**(x; f, p)*Å

*c**n,0**(x)*
*c**n,1**(x)*

:
*c**n,p01**(x)*

,

* p(x)*Å

*Å*

**p(x; p)***x*

^{0}

*x*

^{1}:

*x*

^{p01}.

**It is easy to show that the matrix M(f, p) is nonsingular for all choices of f, p.**

**From definition (2.5), it follows immediately that M**Å **B**^{r}**V, where B is the lower**
triangular matrix with entries

*b**is* Å

### S

^{s}^{i}### D

^{M}

^{i0s}^{,}

^{0}

^{£}

^{s}^{£}

^{i,}**and V is the Vandermonde matrix with entries**^{£}*sj* Å*j** ^{s}*. Hence,

**FIG. 1.** *Construction of b(x; p, x*0).

**det(M)**Å**det(B)**^{r}**det(V)**Å*M** ^{p}*0(

### ∏

^{n01}*kÅ1*

*k!)* x0.

From

**c***n**(x)*Å **M**^{01}**H**^{01}*n* * p(x)* Å

*h*

^{01/2}

**M**

^{01}

**p(x/h),***we observe that each c**n,j**(x) is a polynomial of degree at most (p* 0*1) in (x/h), and*
that

*c**n,j**(x)* Å*h*^{01/2}*c**0,j**(x/h)*Å2^{n/2}*c**0,j*(2^{n}*x).* (3.9)
*To put this approach into a wavelet-like setting, we select a unit interval I* Å*[x*0,
*x*0/*1) for some (as yet undetermined) point x*0*. Scaled and shifted versions of I are*
*denoted by I**n,s*Å*[(x*0 /*s)h, (x*0/*s*/*1)h).*

*We restrict the use of formula (3.7) to the interval I**n,0**. Values of f on a shifted*
*interval I**n,s**are recovered by using the same coefficients c**n,j*on a shifted set of scaling
function coefficients,

*f (x)*É^{p01}

### ∑

*jÅ0*

*f**n,j/s**c**n,j**(x*0*sh)*Å^{s/p01}

### ∑

*jÅs*

*f**n,j**c**n,j0s**(x*0*sh)* *for x*√*I**n,s*. (3.10)

We can combine (3.9) and (3.10) into the desired form (3.2), (3.3) by defining

*b(x; p, x*0)Å

*c**0,p01**(x* /*p*0*1; p)* *if x*√*[x*00*p*/*1, x*00*p*/2),
:

*c*0,1*(x*/*1; p)* *if x*√*[x*00*1, x*0),
*c*0,0*(x; p)* *if x*√*[x*0*, x*0/1).

(3.11)

*By construction, each b(x; p, x*0*) is a piecewise polynomial of degree p* 0 1 with
*support [x*0 0*p*/*1, x*0 /1) (see Fig. 1).

*Formulas for reconstructing the derivatives of f can be obtained by an analogous*

construction. We use the notation b* ^{[r]}* for the function used to reconstruct values of

*the r th derivative; f*

^{(r)}*denotes the rth derivative of f.*

*We again start by looking for coefficients c*^{[r]}*n,j*so that

*f*^{(r)}*(x)*Å^{p01}

### ∑

*jÅ0*

*f**n,j**c*^{[r]}*n,j**(x)* *for f (x)* Å*x*^{k}*, k*Å*0, . . . , p* 01.

This leads to the matrix equation

**H***n***Mc***n*^{[r]}*(x)*Å**p**^{(r)}*(x),*
from which we easily obtain

**c**^{[r]}*n**(x)* Å*h*^{0r}**c**^{(r)}*n**(x)*
and, therefore,

b^{[r]}*(x; p, x*0)Å

*c*^{(r)}*0,p01**(x* /*p*0*1; p)* *if x*√*[x*00*p*/*1, x*00*p*/2),
:

*c** ^{(r)}*0,1

*(x*/

*1; p)*

*if x*√

*[x*00

*1, x*0),

*c*

*0,0*

^{(r)}*(x; p)*

*if x*√

*[x*0

*, x*0/1).

(3.12)

Thus,

b^{[r]}*(x; p, x*0)Åb^{(r)}*(x; p, x*0)

in the sense of b^{[r]}*being the rth derivative from the right of b, i.e., ignoring possible*
*discontinuities at the knots (x*0/ *j), j* √ Z. b* ^{[r]}* is therefore a piecewise polynomial

*of degree (p*0

*r*01).

*Remarks.*

1. It is well known that Vandermonde matrices have condition numbers which increase rapidly with the size of the matrix. We were not able to avoid this phenomenon by alternative approaches to the derivation of b.

*However, b functions for larger p do not appear to offer any real advantages over*
*lower order ones. The lower order functions (p* Å 3 or 4) already give excellent
reconstruction results. Some of the higher order functions are oscillatory and thus
unsuitable, because of cancellation problems during the reconstruction. The others
*look very similar to b functions for smaller p (see Fig. 4), but require a wider support.*

*For small values of p, we observed no difficulties in the numerical calculation of*
*the coefficients. Even for p*Å*10, we obtained the coefficients of the polynomials c**n,j*

to 10 decimals of accuracy, using the Bjo¨rck – Pereyra algorithm [1] in 15-decimal arithmetic.

2. If the dual scaling function f_{H} is a spline, does b equal f_{H}? The answer is: possibly
yes, but usually not.

*Assume that for some integer p, f*_{H} *is a piecewise polynomial of degree p*01, has
*support of length p and can be used to reconstruct all polynomials up to order p*0
1. These conditions are satisfied if f_{H} *is the B-spline of order p, for example.*

*Then, if we use the same p, and if we choose x*0 so that the supports of b and f_{H}
match, the polynomial pieces of both b and f_{H} satisfy (3.8) and by uniqueness f_{H} must
*equal b. For any other choice of p or x*0, b will be different from f_{H}.

Figure 3 shows the dual scaling function f_{H} *and b(x; p) for p* Å 2, 3, 4, for the
*Daubechies – Cohen spline wavelet with (N, N˜ )*Å*(3, 3). b(x; 3) is equal to f*_{H} if we
*choose x*0 Å *1, which makes b continuous. b(x; 2) and b(x; 4) are different from*
f_{H}*, and so is b(x; 3) for other choices of x*0.

3. A natural question to ask is whether b or b* ^{[r]}*satisfy recursion relations analogous
to (2.1). This is not the case in general, as can be verified from the examples in
Section 6.

4. We have not used any special properties of f other than*M*0Å

### *

^{f(x)dx}^{x}

^{0.}

Our construction is valid in any setting where data of the form^{»}*f, f**n,j*…are given.

*5. We will discuss the choice of x*0below. At this point we would just like to point
*out that it is not required to choose the same x*0for b and b^{[r]}*. If x*0*and x** ^{[r]}*0 are chosen
differently, b

*will be different from b*

^{[r]}*(see Figs. 2c, d).*

^{(r)}**4. ERROR ESTIMATES**

*Fix p and x*0*. For x*√*I* Å *[x*0*, x*0/ 1), the pointwise reconstruction error for the
*monomial x*^{k}*at resolution h* Å2^{0}is

*e**k**(x)* Å*x** ^{k}*0

^{p01}### ∑

*jÅ0*

»*x** ^{k}*, f

*0,j*…

*c*

*0,j*

*(x)*Å

*x*

*0*

^{k}

^{p01}### ∑

*jÅ0*

*M**k,j**c**0,j**(x).*

*By (3.8), e**k*å *0 for k*Å*0, . . . , p*0*1. For k*§*p, e**k* *is a polynomial of degree k.*

*If f*√*C*^{p}*, then for x*√*I**n,0*,

*f (x)*Å ^{p01}

### ∑

*jÅ0*

*f** ^{(j)}*(0)

*j!* *x** ^{j}*/

*f*

*(j)*

^{(p)}*p!* *x** ^{p}* for some j√

*I*

*n,0*.

*The reconstruction error for f at the point x is*

*f (x)*0^{p01}

### ∑

*jÅ0*

»*f, f**n,j*…*c**n,j**(x)*Å^{p01}

### ∑

*jÅ0*

*f** ^{(j)}*(0)

*j!* *h*^{j}*e**j**(x/h)* /*f** ^{(p)}*(j)

*p!* *h*^{p}*e**p**(x/h) for some j*√*I**n,0*

Å*O(h** ^{p}*), (4.1)

*since f*^{(p)}*and e**p**are continuous in I**n,0*.

*For x*√*I**n,s*we obtain a corresponding result using Taylor series expansion around
*the point sh.*

*If f*√*C*^{p/1}*, we can take the Taylor series expansion one step further: For x*√*I**n,s*,

*f (x)*0^{p01}

### ∑

*jÅ0*

*f**n,j**c**n,j**(x)*Å *f*^{(p)}*(sh)*

*p!* *h*^{p}*e**p**((x/h)*0*s)* /*O(h** ^{p/1}*). (4.2)

*If s is a zero of e*

*p*

*inside I, the local error at the scaled and translated points (s*/

*s)h is O(h*

*).*

^{p/1}*Remarks.*

*1. It is conceivable that the term O(h** ^{p/1}*) in (4.2) could also vanish at the same
point s, but that requires specific relationships between the higher moments of f.

Unless we are in a position to prescribe the moments of f, we would not expect that to happen.

2. Let s*k**, k*Å*1, 2, . . . , denote the real zeros of e**p**. At these points, we potentially*
get higher order accuracy, but only if s*k**lies inside I*Å*[x*0*, x*0/1). See the numerical
examples in Section 6 for illustration.

*3. It is possible to estimate the error at the point x* √*I**n,s* based on Taylor series
*expansions around a point different from sh, for example around x itself. Such attempts*
lead back to the same results presented here, after laborious calculation.

Differentiation commutes with all the operations leading to the error estimates (4.1)
*and (4.2). In particular, if f*√*C*^{p}*and r*õ *p, then for x*√*I**n,0*

*f*^{(r)}*(x)* Å^{p0r01}

### ∑

*jÅ0*

*f** ^{(j/r)}*(0)

*j!* *x** ^{j}*/

*f*

*(j)*

^{(p)}*(p* 0*r)!x** ^{p0r}* for some j√

*I*

*n,0*

Å^{p01}

### ∑

*jÅ0*

*f** ^{(j)}*(0)

*j!*

*d*^{r}

*dx*^{r}*x** ^{j}*/

*O(h*

*),*

^{p0r}*f*

*n,j*Å

^{p01}### ∑

*iÅ0*

*f** ^{(i)}*(0)

*i!* *h*^{i/(1/2)}*M**i,j*/*O(h** ^{p/(1/2)}*) as before,

*c*

^{[r]}*n,j*Å

*h*

^{0r}*d*

^{r}*dx*^{r}*c**n,j**(x),*

*which leads to an error estimate at x*√*I**n,s*of
*f*^{(r)}*(x)*0^{p01}

### ∑

*jÅ0*

*f**n,j**c*^{[r]}*n,j**(x)* Å^{p01}

### ∑

*jÅ0*

*f*^{(j)}*(sh)*
*j!*

*d*^{r}

*dx*^{r}*[h*^{j}*e**j**((x/h)*0*s)]*

/*f*^{(p)}*(sh)*
*p!*

*d*^{r}

*dx*^{r}*[h*^{p}*e**p**((x/h)* 0*s)]* /*O(h** ^{p0r/1}*)

Å*O(h** ^{p0r}*). (4.3)

Again, the order of convergence is one level higher at certain points. These points
s^{[r]}*k* *are the zeros of e*^{[r]}*p* *in I**n,s*and are, in general, different from s*k*.

**5. SMOOTHNESS PROPERTIES OF THE b FUNCTIONS**

*In Sections 3 and 4, we had to perform the calculations for an arbitrary level n*
*in order to obtain the correct powers of h. Since smoothness is independent of*

*scaling, it suffices to take n*Å *0, h*Å 2^{0}in this section, and we use the notation
*c**j**instead of c**0,j*.

T^{HEOREM}5.1. *The function b*^{[r]}*(x; x*0*) is continuous everywhere if it is continuous*
*at one of the endpoints of its support.*

*Proof.* *We will write out the proof for r* Å 0. As in the derivation of (4.3),
the general proof is then obtained by differentiating everything, since differentiation
commutes with all operations we use.

Continuity of b means (see (3.11) and Fig. 1)
0Å*c*0*(x*0/1),

*c**j**(x*0) Å*c**j/1**(x*0/1), *j*Å*0, . . . , p*02,
*c**p01**(x*0) Å0.

*Without loss of generality, assume continuity at the left endpoint, that is, c**p01**(x*0) Å
**0. Since M is nonsingular, it suffices to prove**

**M**
0
*c*0*(x*0)

:
*c**p02**(x*0)

Å**M**

*c*0*(x*0/1)
*c*1*(x*0/1)

:
*c**p01**(x*0/ 1)

.

*The kth entry on the right-hand side is (x*0 /1)^{k}*by (3.8). Using (2.6), the kth entry*
on the left is

*p02*

### ∑

*jÅ0*

*M**k,j/1**c**j**(x*0) Å^{p01}

### ∑

*jÅ0*

*M**k,j/1**c**j**(x*0) *(since c**p01**(x*0)Å0)

Å^{p01}

### ∑

*jÅ0*

### H

*lÅ0*

^{∑}

^{k}### S

^{k}^{l}### D

^{M}

^{l,j}### J

^{c}

^{j}

^{(x}^{0}

^{)}

^{Å}

*lÅ0*

^{∑}

^{k}### S

^{k}^{l}### D

^{x}

^{l}^{0}

^{Å}

^{(x}^{0}

^{/}

^{1)}

^{k}^{.}

^{j}

*For even p, c**p01*is a polynomial of odd degree, so there always exists a choice of
*x*0which makes b continuous. It is not obvious whether that is always possible for
*odd p, but we were able to find such x*0in all the examples we considered.

*As illustrated in the numerical examples below, x** ^{[r]}*0 should be chosen so that the

*leading error term e*

^{[r]}*p*

*is small over the interval [x*

*0*

^{[r]}*, x*

*0 / 1), whether or not this makes b*

^{[r]}

^{[r]}*continuous. However, in all examples we tried we were able to find x*

*0*

^{[r]}which produces a continuous b* ^{[r]}* and small error simultaneously.

The following lemma, together with Theorem 5.1, states that the starting points
*x*0* ^{[r]}*which make b

^{[r]}*(x; p*/1) continuous are the same as the points s

*of (potentially) higher accuracy for b*

^{[r]}

^{[r]}*(x; p).*

LEMMA5.2. *For any r*õ*p, if c*^{[r]}*p**(s; p*/1) Å*0, then*

*c*^{[r]}*j* *(s; p*/1)Å*c*^{[r]}*j* *(s; p),* *j* Å*0, . . . , p*01, (5.1)

*and*

*e*^{[r]}*p*(s)Å *d*^{r}

*ds** ^{r}*s

*0*

^{p}

^{p01}### ∑

*jÅ0*

*M**p,j**c*^{[r]}*j* *(s; p)*Å0. (5.2)
*Proof.* *We will show the proof for r* Å*0 only; the proof for general r proceeds*
along identical lines, with derivatives inserted all over the place.

*Assume c**p**(s; p*/1)Å0, then

* M(p*/

*/1) Å*

**1)c(s; p***M**0,p*

* M(p)* :

*M**p01,p*

*M**p,0* ??? *M**p,p01* *M**p,p*

*c*0*(s; p*/1)
:
*c**p01**(s; p*/1)

0

(5.3)

Å
s^{0}

:
s^{p01}

s* ^{p}*
.

Thus,

**M(p)**

## S

^{c}^{c}

^{p01}^{0}

^{(s; p}^{(s; p}^{:}

^{/}

^{/}

^{1)}

^{1)}

## D

^{Å}

## S

^{s}

^{s}

^{:}

^{p01}^{0}

## D

which implies (5.1) by uniqueness of the solution of (3.8).

Equation (5.2) is the last equation in (5.3), using (5.1). j

We summarize the results of the last three sections in the following theorem.

THEOREM 5.3. *For any choice of f, p, r*õ *p, x** ^{[r]}*0

*, there exists a function b*

^{[r]}*(x)*Åb

^{[r]}*(x; f, p, x*

*0*

^{[r]}*) so that for f*√

*C*

^{p}*and for any x*√R

*,*

É*f*^{(r)}*(x)* 0*B*^{[r]}*n* *f (x)*^{É} Å*O(h** ^{p0r}*),

*where*

*h*Å2* ^{0n}*,

*B*

^{[r]}*n*

*f*Å

### ∑

^{`}

*jÅ0*`

*f**n,j*b^{[r]}*n,j*,
b^{[r]}*n,j**(x)* Å2* ^{n/2}*b

*(2*

^{[r]}

^{n}*x*0

*j).*

b^{[r]}*(x; f, p, x** ^{[r]}*0

*) has support [x*

*0 0*

^{[r]}*p*/

*1, x*

*0 /*

^{[r]}*1] and is a polynomial of degree at*

*most p*0

*r*0

*1 on each subinterval [x*

*0 /*

^{[r]}*j, x*

*0 /*

^{[r]}*j*/

*1], j*√Z

*.*

**TABLE 1**

**Choices of Starting Points x**^{[r]}**0** **Which Make b**^{[r]}**Continuous and Points s**^{[r]}**of Potentially**
**Higher Accuracy for (a) Orthogonal Daubechies Wavelet with Two Vanishing Moments**
**(b) Biorthogonal Daubechies – Cohen Spline Wavelet with (N, N****˜ ) Å (3, 3)**

(a) (b)

*x** ^{[r]}*0 s

^{[r]}*x*

*0 s*

^{[r]}

^{[r]}*p*Å2 *r*Å0 (30 1/2 1

q

3)/2 (30

q

3)/2 (50

q

3)/2

*r*Å1 (40 1

q

3)/2

*p*Å3 *r*Å0 (30 0.57706 1 1

q

3)/2

**FIG. 2.** Examples of b* ^{[r]}*for Daubechies wavelet with two vanishing moments: (a) scaling function

*f; (b) b (x; p*Å

*2, x*0Å

*M*1

*); (c) b (x; p*Å

*3, x*0Å

*M*1/1); (d) b

^{[1]}

*(x; p*Å

*3, x*0Å

*M*1/

^{1}

_{2}

*); (e) b(x; p*Å

*3, x*0Å0).

*e*2(s)Å s^{2}0*M*2,0*c*0(s)0 *M*2,1*c*1(s)Å0,
which leads to

sÅ^{1}_{2}(2*M*1 /1{

q

4*M*204*M*^{2}1/1).

If*M*2Å *M*^{2}1, which happens for all orthogonal wavelets with at least two vanishing
moments, then sÅ*M*1or*M*1/1.

b^{[1]} is always a step function and is never continuous. The leading error term for
*B*^{[1]}*n* vanishes if

**FIG. 2—Continued**

**FIG. 3.** Examples of b^{[r]}*for biorthogonal Cohen – Daubechies spline wavelet with (N, N)*Å(3, 3): (a)
dual scaling function f_{H}*; (b) b(x; p*Å*2, x*0Å^{1}_{2}*); (c) b(x; p*Å*3, x*0Å*1); (d) b(x; p*Å*4, x*0Å^{3}_{2}).

*e**2(s^{[1]}) Å2s^{[1]}0 *M*2,0*c**0(s^{[1]})0 *M*2,1*c**1(s^{[1]})Å 0,
which leads to

s^{[1]} Å*M*1/^{1}_{2}.
*For p*Å3, b is continuous if

*x*0Å^{1}_{2}(2*M*1/1{

q

4*M*204*M*^{2}1/ 1),
and b^{[1]}is continuous if

*x*^{[1]}0 Å*M*1/^{1}_{2}.

**FIG. 3—Continued**

These are the s, s^{[1]} *for p*Å 2 (see remark after Lemma 5.2). b^{[2]} is always a step
function and is never continuous.

*The leading error term for B**n*vanishes if

*e*3(s)Ås^{3}0*M*3,0*c*0(s)0*M*3,1*c*1(s)0*M*3,2*c*2(s)Å0.

The resulting cubic equation cannot be solved in closed form in general. If *M*2 Å

*M*^{2}1,*M*3Å*M*^{3}1, which happens in the case of higher order coiflets, then sÅ*M*1,*M*1

/1 or*M*1/2.

*The leading error term for B*^{[1]}*n* vanishes for

s^{[1]} Å1/*M*1{^{q}1/30*M*^{2}1/*M*2.
*The leading error term for B*^{[2]}*n* vanishes for

s^{[2]} Å1/*M*1.

**FIG. 4.** *b for large p, for Daubechies wavelet with two vanishing moments: (a) b(x; p*Å*5); (b) b(x;*

*p*Å*10); (c) b(x; p*Å15).

**TABLE 2**

**Comparison of Errors at Several Points x and Several Resolutions h, for Reconstruction****of sin x Using b(x; p Å 3, x****0** **Å 1.633975), Based on Daubechies Wavelet with Two**
**Vanishing Moments**

*x* *h*Å2^{0} *h*Å2^{01} *h*Å2^{02}

0.25 0.0020729566 0.0005046129 0.0000621160

. . .

0.50 0.0040526588 0.0004775991 0.0000575524

. . .

1.00 0.0031835356 0.0003380129 0.0000381963

*Note. The underlined values show the predicted convergence behavior O(h*^{3}).

As particular examples, we consider two types of wavelets:

• The orthogonal Daubechies wavelet with 2 vanishing moments [3, p. 194 ff]. Its first few moments are

*M*0Å1, *M*1Å(30

q

3)/2 É0.633975, *M*2Å3(20

q

3)/2 É0.40192.

• *The biorthogonal Daubechies – Cohen spline wavelet with (N, N˜ )* Å (3, 3) (see
[2, p. 543] or [3, p. 271 ff]). Its first few moments are

*M*0Å1, *M*1Å^{1}_{2}, *M*2Å0, *M*3Å 0^{1}_{4},

*The points x** ^{[r]}*0 and s

^{[r]}*for p*Å2 and 3 are summarized in Table 1. The dual scaling functions and some of the functions b, b

^{[1]}are shown in Figs. 2, 3, and 4.

*6.2. Reconstruction Examples*

The numerical examples were done on a DECstation 5000/133, using Matlab 4.2c.

*The scaling function coefficients f**n,j* were calculated numerically to accuracy 10^{012},
all subsequent calculations were done by Matlab in double precision.

*Most of the examples are based on the reconstruction of the function f (x)*Å*sin x*
on the interval [0, 3].

First, we illustrate the error estimates (4.1) through (4.3), using the Daubechies
*wavelet with two vanishing moments and accuracy level p*Å3. The underlined values
*in Table 2 are approximately proportional to f*-*(0)h*^{3} and decrease by the expected
*factor of 8 at each level; the values in each row (errors at a fixed point x for different*
resolutions) are not directly comparable.

**TABLE 3**

**Comparison of Maximum Errors at Various Resolutions**
**for Reconstruction of sin x or Its Derivative on [0, 3]**

Maximum errors Reconstruction Errors

using evaluated for *h*Å2^{0} *h*Å2^{01} *h*Å2^{02} *h*Å2^{03}

(a) f *Dx*Å1/64 0.2230931 0.0597919 0.0154932 0.0038994

b *Dx*Å1/64 0.0418520 0.0054484 0.0006783 0.0000842

s1Å2.565179 0.0179483 0.0012485 0.0000801 0.0000050
s2Å1.759679 0.0072115 0.0004933 0.0000316 0.0000020
s3Å0.577066 0.0162341 0.0012062 0.0000898 0.0000079
b^{[1]} *Dx*Å1/64 0.1420792 0.0402152 0.0103265 0.0025981

(b) f *Dx*Å1/64 0.2230931 0.0597919 0.0154932 0.0038994

b *Dx*Å1/64 0.3933545 0.0531412 0.0067300 0.0008467

(c) f *Dx*Å1/64 0.0983246 0.0084478 0.0005585 0.0000354

b *Dx*Å1/64 0.0993215 0.0139191 0.0017952 0.0002273

Ratios

(a) f *Dx*Å1/64 3.7311581 3.8592332 3.9732578

b *Dx*Å1/64 7.6815296 8.0325838 8.0536774

s1Å2.565179 14.3753306 15.5893594 15.9203390 s2Å1.759679 14.6201103 15.6096820 15.8281089 s3Å0.577066 13.4590786 13.4291017 11.4237320

b^{[1]} *Dx*Å1/64 3.5329707 3.8943582 3.9746441

(b) f *Dx*Å1/64 3.7311581 3.8592332 3.9732578

b *Dx*Å1/64 7.4020651 7.8961797 7.9484083

(c) f *Dx*Å1/64 11.6390266 15.1266675 15.7609234

b *Dx*Å1/64 7.1356205 7.7535307 7.8987923

*Note. The errors were calculated either at equally spaced points with spacing Dx, or at all scaled and*
translated versions of the points s*i*. The second part of the table shows ratios between adjacent columns
*in the first part. (a) Using continuous b(x; p*Å*3, x*0Å1.633975) and b^{[1]}*(x; p*Å*3, x*0Å1.133975) based
*on Daubechies wavelet with two vanishing moments. (b) Using discontinuous b(x; p*Å*3, x*0Å0) based
on Daubechies wavelet with two vanishing moments. (c) Using continuous b^{[1]}*(x; p*Å*3, x*0Å1.133975)
*based on the biorthogonal Daubechies-Cohen spline wavelet with (N, N**˜ )*Å(3, 3).

Table 3(a) shows the maximum errors at equally spaced points across the interval
[0, 3] for the original scaling function series (1.1) and for continuous b and b^{[1]}with
*p*Å3. Also shown are the maximum errors at scaled translates of s*k*.

As expected, for general points the ratio of errors at different resolutions approaches
4 for f and b^{[1]}*and 8 for b. This indicates errors of order O(h*^{2}*) and O(h*^{3}), respectively.

For the points s1and s2*, which are inside I, convergence is of order O(h*^{4}), s3lies
*outside I, and the ratio approaches the value of 8, corresponding to the error O(h*^{3})
valid for general points.

*Figure 5 shows the reconstructions themselves, at resolution h* Å 2^{0}. Figure 6

**FIG. 5.** *Reconstruction of sin x at resolution h*Å2^{0}: Dotted line, original function; dashed line, series
*using Daubechies wavelet with two vanishing moments; solid line, series using b(x; p*Å*3, x*0Å*M*1/1).

compares the reconstruction errors for f and for b, as well as the error at scaled
translates of s1*at resolution h* Å2^{04}; note the logarithmic scale on the vertical axis.

*Second, we illustrate the effect of the choice of x*0. Part (b) of Table 3 is comparable
*to the first two lines of part a, except that we have chosen x*0Å0. For this choice, b
is not continuous (see Fig. 2e). The errors still decay at the expected rate, but they
are now much larger.

*To understand the influence of x*0 *better, consider the graph of e*3 in Fig. 7. The
*error depends on the size of e*3*over the interval I*Å*[x*0*, x*0/1). In general, we would

**FIG. 6.** *Pointwise reconstruction error in reconstruction of sin x at resolution h*Å2^{04}: Dashed line,
*error using Daubechies wavelet with two vanishing moments; solid line, error using b(x; p*Å*3, x*0Å*M*1

/1); circled points, error at scaled translates of s1É2.56518.

**FIG. 7.** *Leading error term e*3*(x) for b(x; p*Å3), based on Daubechies wavelet with two vanishing
moments. Zeros of the function (i.e., points s*j*points) are circled.

*expect e*3*to be small if I is located somewhere in the region between the zeros of e*3

*(such as x*0Å*1.633975) and larger if I is farther away from the region of zeros (such*
*as x*0Å0).

The size of the error does not depend on whether b is continuous or not, but it is
*esthetically more pleasing to have a continuous approximating function if an x*0can
*be found in the region where e**p*is small.

Third, we demonstrate that our results work equally well for biorthogonal wavelets.

Part (c) of Table 3 shows results comparable to the first two lines of part (a), for the

**FIG. 8.** *Reconstruction of sin x*/0.51*sin(2.5(x*01))/0.71*sin(6(x*/*1)) at resolution h*Å2^{0}:
Dotted line, original function; dashed line, series using Daubechies wavelet with two vanishing moments;

*solid line, series using b(x; p*Å*3, x*0Å*M*1/1).

**FIG. 9.** *Pointwise reconstruction error in reconstruction of sin x* /0.51 *sin(2.5(x*01))/ 0.71
*sin(6(x*/*1)) at resolution h*Å2^{0}*: Solid line, error using b(x; p*Å*3, x*0Å*M*1/1); circled points, error
at scaled translates of s1É2.56518.

*biorthogonal Daubechies – Cohen spline wavelet with (N, N˜ )*Å(4, 2) (see [3, p. 271
ff]). In this case, our reconstruction is less accurate, but faster, than the scaling function
series.

Experiments with other wavelets and with other smooth and slowly varying func- tions yield similar results to the ones presented so far.

Finally, Figs. 8 and 9 illustrate what happens in the case of a highly oscillating function; the results for functions with singularities are similar.

*As shown in Fig. 8, the reconstruction B**n**f is smoother looking, but not any closer*
*to the original curve than the scaling function series P˜*_{n}*f. However, the accuracy of*
*B**n**f is not any worse than the accuracy of P˜*_{n}*f, either.*

Figure 9 illustrates that the error at the scaled translates of s*k* is not noticeably
smaller than at other points. Obviously, the error in the case of a highly oscillating
function is not described all that well by the leading term.

**7. SUMMARY**

*For situations where it is desired to recover point values of a function f to high*
*accuracy from its scaling function coefficients f**n,j*, we propose to replace the standard
scaling function series

*P*H^{n}*f* Å

### ∑

^{`}

*jÅ0`*

»*f, f**n,j*…f_{H}*n,j*

by a corresponding series

*B**n**f*Å

### ∑

^{`}

*jÅ0*`

»*f, f**n,j*…b*n,j*,

or more generally,

*f*^{(r)}*(x)* É*B*^{[r]}*n* *f*Å

### ∑

^{`}

*jÅ0`*

»*f, f**n,j*…b^{[r]}*n,j*.

*Theoretical estimates and numerical experiments indicate that the new series B**n**f has*
*accuracy at least as good as the scaling function series for all functions f and performs*
*significantly better for smooth, slowly varying f.*

*In addition, derivatives of f can be recovered by using the related series B*^{[r]}*n**f for*
which there is in general no counterpart using derivatives of scaling functions.

**ACKNOWLEDGMENT**

The authors thank the (anonymous) referee for many helpful hints which led to improved clarity and simplified proofs.

**REFERENCES**

1. A* ˚ . Bjo¨rck and V. Pereyra, Solution of Vandermonde systems of equations, Math. Comp. 24 (1970),*
893 – 903.

2. A. Cohen, I. Daubechies, and J.-C. Feauveau, Biorthogonal bases of compactly supported wavelets,
**Comm. Pure Appl. Math. 45 (1992), 485 – 560.**

3. I. Daubechies, ‘‘Ten Lectures on Wavelets,’’ CBMS-NSF Regional Conference Series in Applied Mathematics, Vol. 61 SIAM, Philadelphia, 1992.

*4. S. Kelly, M. Kon, and L. Raphael, Pointwise convergence of wavelet expansions, Bull. Amer. Math.*

**Soc. 30 (1994), 87 – 94.**

*5. S. G. Mallat, Multiresolution approximations and wavelet orthonormal bases of L*^{2}(_{R}*), Trans. Amer.*

**Math. Soc. 315 (1989), 69 – 87.**

6. W. Sweldens and R. Piessens, Quadrature formulae and asymptotic error expansions for wavelet
**approximations of smooth functions, SIAM J. Numer. Anal. 31 (1994), 1240 – 1264.**

7. G. G. Walter, ‘‘Wavelets and Other Orthogonal Systems with Applications,’’ CRC Press, Boca Raton, FL, 1994.