See discussions, stats, and author profiles for this publication at: [Link]

net/publication/2580451

Wavelets

Article in Acta Numerica · September 1998

DOI: 10.1017/S0962492900002233 · Source: CiteSeer

CITATIONS READS

66 172

2 authors:

Ronald DeVore Bradley J. Lucier

Texas A&M University Purdue University West Lafayette
274 PUBLICATIONS 24,575 CITATIONS 67 PUBLICATIONS 3,328 CITATIONS

SEE PROFILE SEE PROFILE

All content following this page was uploaded by Bradley J. Lucier on 07 December 2013.

The user has requested enhancement of the downloaded file.

Ada Numerica (1991), pp. 1-56

Wavelets*
Ronald A. DeVore
Department of Mathematics
University of South Carolina,
Columbia, SC 29208 USA
E-maU: devore@[Link]
Bradley J. Lucier
Department of Mathematics
Purdue University,
West Lafayette, IN 47907 USA
E-mail: lucier@[Link]

CONTENTS
1 Introduction 1
2 The Haar wavelets 4
3 The construction of wavelets 13
4 Fast wavelet transforms 37
5 Smoothness spaces and wavelet coefficients 40
6 Applications 44
References 54

1. Introduction
The subject of 'wavelets' is expanding at such a tremendous rate that it is
impossible to give, within these few pages, a complete introduction to all
aspects of its theory. We hope, however, to allow the reader to become suf-
ficiently acquainted with the subject to understand, in part, the enthusiasm
of its proponents toward its potential application to various numerical prob-
lems. Furthermore, we hope that our exposition can guide the reader who
wishes to make more serious excursions into the subject. Our viewpoint is
biased by our experience in approximation theory and data compression; we
warn the reader that there are other viewpoints that are either not repre-
sented here or discussed only briefly. For example, orthogonal wavelets were
developed primarily in the context of signal processing, an application upon

* This work was supported in part by the National Science Foundation (grants DMS-
8922154 and DMS-9006219), the Air Force Office of Scientific Research (contract 89-
0455-DEF), the Office of Naval Research (contracts N00014-90-1343, N00014-91-J-1152,
and N00014-91-J-1076), the Defense Advanced Research Projects Agency (AFOSR con-
tract 90-0323), and the Army High Performance Computing Research Center.
 R. A. D E V O R E AND B. J. LUCIER

4>

-1 0 j - 1 j +1
_ 2*
2* 2*
Fig. 1. An example of functions <f> and ^(2*- — j).

which we touch only indirectly. However, there are several good expositions
(e.g. Daubechies (1990) and Rioul and Vetterli (1991)) of this application.
A discussion of wavelet decompositions in the context of Littlewood-Paley
theory can be found in the monograph of Frazier et al. (1991). We shall also
not attempt to give a complete discussion of the history of wavelets. Histor-
ical accounts can be found in the book of Meyer (1990) and the introduction
of the article of Daubechies (1990). We shall try to give sufficient historical
commentary in the course of our presentation to provide some feeling for
the subject's development.
The term 'wavelet' (originally called wavelet of constant shape) was intro-
duced by J. Morlet. It denotes a uni-variate function ty (multi-variate wave-
lets exist as well and will be discussed subsequently), defined on R, which,
when subjected to the fundamental operations of shifts (i.e. translation by
integers) and dyadic dilation, yields an orthogonal basis of L2(R). That is,
the functions V'[Link] := 2*/2V>(2*- — j), j , f c 6 Z , form a complete orthonormal
system for L2OR). In this work, we shall call such a function an orthogo-
nal wavelet, since there are many generalizations of wavelets that drop the
requirement of orthogonality. The Haar function H := X[o,i/2) ~ X[i/2,i)>
which will be discussed in more detail in the section that follows, is the sim-
plest example of an orthogonal wavelet. Orthogonal wavelets with higher
smoothness (and even compact support) can also be constructed. But before
considering that and other questions, we wish first to motivate the desire
for such wavelets.
We view a wavelet ip as a 'bump' (and think of it as having compact
support, though it need not). Dilation squeezes or expands the bump and
translation shifts it (see Figure 1). Thus, tpj^ is a scaled version of xp centred
at the dyadic integer j2~k. If k is large positive, then tpj^ is a bump with
small support; if k is large negative, the support of ipj^ is large.
 WAVELETS 3

The requirement that the set {V"j,fc}j,fcez forms an orthonormal system

means that any function / G L2O&) can be represented as a series
/ =

with {/, g) := JRfgdx the usual inner product of two L 2 W functions. We

view (1.1) as building up the function / from the bumps ipj,k- Bumps
corresponding to small values of k contribute to the broad resolution of / ;
those corresponding to large values of k give finer detail.
The decomposition (1.1) is analogous to the Fourier decomposition of a
function / G Z<2(T) in terms of the exponential functions e^ :=eik', but there
are important differences. The exponential functions e*; have global support.
Thus, all terms in the Fourier decomposition contribute to the value of /
at a point x. On the other hand, wavelets are usually either of compact
support or fall off exponentially at infinity. Thus, only the terms in (1.1)
corresponding to Vj,fc with j2~k near x make a large contribution at x. The
representation (1.1) is in this sense local. Of course, exponential functions
have other important properties; for example, they are eigenfunctions for
differentiation. Many wavelets have a corresponding property captured in
the 'refinement equation' for the function <j> from which the wavelet ijj is
derived, as discussed in Section 3.1.
Another important property of wavelet decompositions not present di-
rectly in the Fourier decomposition is that the coefficients in wavelet de-
compositions usually encode all information needed to tell whether / is in a
smoothness space, such as the Sobolev and Besov spaces. For example, if ip
is smooth enough, then a function / is in the Lipschitz space Lip(o!, Loc(R)),
0 < a < 1, if and only if

sup2* (a+ 2 ) |(/,t/' : , i fc)| (1.2)

is finite, and (1.2) is an equivalent semi-norm for Lip(a, Loo(R)).

All this would be of little more than theoretical interest if it were not for
the fact that one can efficiently compute wavelet coefficients and reconstruct
functions from these coefficients. Such algorithms, known as 'fast wavelet
transforms' are the analogue of the Fast Fourier Transform and follow simply
from the refinement equation mentioned earlier.
In many numerical applications, the orthogonality of the translated di-
lates tpjtk is not vital. There are many variants of wavelets, such as the
pre-wavelets proposed by Battle (1987) and the ^-transform of Frazier and
Jawerth (1990), that do not require orthogonality. Typically, for a given
function ip, one wants the translated dilates V'[Link]) h & € Z, to form a stable
basis (also called a Riesz basis) for Z,2(K). This means that each / € L 2 W
has a unique series decomposition in terms of the [Link], and that the £2
4 R. A. DEVORE AND B. J. LUCIER

norm of the coefficients in this series is equivalent to ||/||L 2 (R) (this will be
discussed in more detail in Section 3.1). In other applications, when ap-
proximating in Li(R), for example, one must abandon the requirement that
V'[Link]) h k € Z, form a stable basis of Z/i(R), because none exists. (The Haar
system is a Schauder basis for Li([0,1]), for example, but the representa-
tion is not Li([0, l])-stable.) For such applications, one can use redundant
representations of / , with ip a box spline, for example.
We have, to this point, restricted our discussion to uni-variate wavelets.
There are several constructions of multi-variate wavelets but the final form
of this theory is yet to be decided. We shall discuss two methods for con-
structing multi-variate wavelets; one is based on tensor products while the
other is truly multi-variate.
The plan of the paper is as follows. Section 2 is meant to introduce the
topic of wavelets by studying the simplest orthogonal wavelets, which are
the Haar functions. We discuss the decomposition of LP(R) using the Haar
expansion, the characterization of certain smoothness spaces in terms of the
coefficients in the Haar expansion, the fast Haar transform, and multi-variate
Haar functions. Section 3 concerns itself with the construction of wavelets.
It begins with a discussion of the properties of shift-invariant spaces, and
then gives an overview of the construction of uni-variate wavelets and pre-
wavelets within the framework of multi-resolution. Later, mention is made
of Daubechies' specific construction of orthonormal wavelets of compact sup-
port. We finish with a discussion of wavelets in several dimensions.
Section 4 examines how to calculate the coefficients of wavelet expan-
sions via the so-called Fast Wavelet Transform. Section 5 is concerned with
the characterization of functions in certain smoothness classes called Besov
spaces in terms of the size of wavelet coefficients. Section 6 turns to nu-
merical applications. We briefly mention some uses of wavelets in nonlinear
approximation, data compression (and, more specifically, image compres-
sion) and numerical methods for partial differential equations.

2. The Haar wavelets

2.1. Overview
The Haar functions are the most elementary wavelets. While they have many
drawbacks, chiefly their lack of smoothness, they still illustrate in the most
direct way some of the main features of wavelet decompositions. For this
reason, we shall consider in some detail their properties which make them
suitable for numerical applications. We hope that the detail we provide at
this stage will render more convincing some of the later statements we make,
without proof, about more general wavelets.
We consider first the uni-variate case. Let H := X[o,i/2) — X[i/2,i) be the
Haar function that takes the value 1 on the left half of [0,1] and the value
 WAVELETS 5

—1 on the right half. By translation and dilation, we form the functions

HiJk := 2k'2H{2k- - j), j , keZ. (2.1.1)
Then, Hj^ is supported on the dyadic interval Ijtt '•= \j2~k, {j + l)2~ fc ).
It is easy to see that these functions form an orthonormal system. In
fact, given two of these functions •#,-,*, #/,*', k' > k and (j, k) ^ (j', k'), we
have two possibilities. The first is that the dyadic intervals /_,-,* and /,•/,# are
disjoint, in which case JR HjtkHj',k' = 0 (because the integrand is identically
zero). The second possibility is that k' > k and Iy^' is contained in one of
the halves J of 7,-^. In this case Hjk is constant on J while Hy^i takes the
values ±1 equally often on its support. Hence, again / R Hj^Hj'tk> — 0.
We want next to show that {•#/,* | j , k € Z} is complete in I ^ W - The
following development gives us a chance to introduce the concept of multi-
resolution, which is the main vehicle for constructing wavelets and which will
be discussed in more detail in the section that follows. Let S := 5 ° denote
the subspace of Z/2(K) that consists of all piecewise-constant functions with
integer breakpoints; i.e. functions in S are constant on each interval \j,j+l),
j € Z. Then S is a shift-invariant space: if S € S, each of its shifts, S(- +k),
k € Z, is also in S. A simple orthonormal basis for <S is given by the shifts
of the function <f> := X[o,i] • Namely, each S € S has a unique representation

S = £ cU)<K • ~ j), (c(j)) e * 2 (z). (2.1.2)

By dilation, we can form a scale of spaces

Sk ~ {S(2k •) | S 6 S}, keZ.
Thus, Sk is the space of piecewise-constant L 2 W functions with breakpoints
at the dyadic integers j2~k. The normalized dyadic shifts <f>jtk'•=2k/2(f>(2k- —
j) = 2k'24>{2k{ • -j2~k)) with step j2~k, j € Z, of the function <f>(2k •), form
an orthonormal basis for Sk. However, to avoid possible confusion, we note
that the totality of all such functions </>jtk is not a basis for the space 2/2 (K)
because there is redundancy. For example, (j> = (<£o,i + 4>i,i)/V2.
Clearly, we have Sk C Sk+1, k G Z, so the spaces <Sfc get 'thicker' as k
gets larger and 'thinner' as k gets smaller. We are interested in the limiting
spaces
S°°:=\Js^ and S"00 := f)Sk, (2.1.3)
since these spaces hold the key to showing that the Haar basis is complete.
We claim that
5°° = L2(R) and 5"°° = {0}. (2.1.4)
The first of these claims is equivalent to the fact that any function in L 2 (K.)
can be approximated arbitrarily well (in the Z/2W norm) by the piecewise-
6 R. A. DEVORE AND B. J. LUCIER

constant functions from Sk provided k is large enough. For example, it is

enough to approximate / by its best L2OR) approximation from Sk. This
best approximation is given by the orthogonal projector Pk from L2W onto
Sk. It is easy to see that

Pkf(x) = j^-j f f, * € Ij,k, j € Z.

\Ij,k\ JIj,k
To verify the second claim in (2.1.4), we suppose that / € f)Sk- Then,
/ is constant on each of (—oo,0) and [0, oo), and since / € L2(R), we must
have / = 0 a.e. on each of these intervals.
Now, consider again the projector Pk from L2OR) onto Sk. By (2.1.4),
Pkf —* / , k —> oo. We also claim that Pkf —• 0, k —* — oo. Indeed, if this
were false, then we could find a C > 0 and a subsequence rrij —• — oo for
which ||Pmj/||L2(R) > C for all m,j. By the weak-* compactness of Z<2(R),
we can also assume that P m j / —* 9, weak-*, for some g € Z«2(K). Now, for
any m € Z, all Pmif arc in Sm for rrij sufficiently large and negative. Since
Sm is weak-* closed, g € <Sm. Hence g e D<Sm implies g = 0 a.e. This gives
a contradiction because by orthogonality

/ |Prai/|2da = / 7PmJdx -» / Jgdx = 0.

J& JtL JU
It follows that each / e L2(R) can be represented by the series
Qk-i:=Pk-Pk-i (2.1.5)
fc€z fcez
because the partial sums, P n / — P-nf, of this series tend t o / as n —• oo.
To complete the construction of the Haar wavelets, we need the following
simple remarks about projections. If Y c X are two closed subspaces of
L2(M) and Px and Py are the orthogonal projectors onto these spaces, then
Q : = Px — Py is t h e orthogonal projector from L2(R) onto X QY, t h e
orthogonal complement of Y in X (this follows from t h e identity PyPx =
P y ) . Thus, t h e operator Qk-i '•= Pk — Pfc-i appearing in (2.1.5) is t h e
orthogonal projector onto Wk~l := Sk Q S*1""1. T h e spaces Wk are t h e
dilates of the wavelet space
W — S^S0. (2.1.6)
k
Since the spaces W , k € Z, are mutually orthogonal, we have W^ -L V0,
j ^ A;, and (2.1.5) shows that L2(M) is the orthogonal direct sum of the Wk:
L2(R) = @Wk. (2.1.7)
fcez
How does the Haar function fit into all this? Well, the main point is that
H and its translates H( • — k) form an orthonormal basis for W. Indeed,
H = 20(2 •) - P o (20(2 •)), which shows that H is in Sl Q S° = W. On the
 WAVELETS 7

other hand, the identities H+</> = 20(2 •) and <p-H = 20(2 • - 1 ) show that
the shifts of <f> together with the shifts of H will generate all the half-shifts
of TJ := <t>(2-). Since the half-shifts of \f2r\ form an orthonormal basis for
S1, the shifts of H must be complete in W.
By dilation, the functions Hj,h, j £ Z, form a complete orthonormal sys-
tem for Wk. Hence, we can represent the orthogonal projector Qk onto Wk
by

jez
Using this in (2.1.5), we have for any / 6 Z<2(R) the decomposition
(2-1-8)
feez jez
In other words, the functions Hjk, j,k € Z, form an orthonormal basis for

2.2. The Haar decomposition in LP(WL)

While the Haar decomposition is initially defined only for functions in L2OR),
it is worth noting that Haar decompositions also hold for other spaces of
functions. In this section, we shall discuss the Haar representation for func-
tions in LP(R), 1 < p < 00. A similar analysis can be given when p = 00
if Loo(R) is replaced by the space of uniformly continuous functions that
vanish at 00, equipped with the Xoo(IR) norm.
If / € LP(R), the Haar coefficients {f,Hjtk) are well defined and we can
ask whether the Haar series (2.1.8) converges in LP(K) to / . To answer
this question, we fix a value of 1 < p < 00 and a k £ Z and examine the
projector Pk, which is initially defined only on L2(K). For any / € 1*2(R), we
have Pkf = Z)/6pfc fiXi where Dk denotes the collection of dyadic intervals
of length 2~k, and where / / := (1/|/|) fjfdx, I € Dk, is the average of
/ over /. In this form, the projector Pk has a natural extension to LP(E)
and takes values in the space Sk(X, LP(R)) of all functions in LP(R) that are
piecewise-constant functions with breakpoints at the dyadic integers j2~k,
j GZ.
In representing Pk on L P (R), it is useful to change our normalization
slightly. We fix a value of p and consider the Lp(R)-normalized characteristic
functions <f>jtk,p := 2k/p<j>(2k--j), (j> := X[o,i], which satisfy / R \(f>j^,p\pdx = 1.
Then,

\ + 1
i€Z
8 R. A. DEVORE AND B. J. LUCIER

Prom Holder's inequality, we find \{f,<t>j,k,j/)\p < //. | / | p d x and so

jez jez
Therefore, Pk is a bounded operator with norm 1 on the space LP(M).
If / e Lp(M), then since Pk is a projector of norm 1,

11/ - JVHMR) = einL IKJ - Pk

W ~ 5 ) U M « ) ^ 2dist(/,5*) Lp(R) . (2.2.1)
It follows that Pkf -> f in L P (R) for each / e L P (R).
On the other hand, consider Pkf as k —> — oo. If / is continuous and
of compact support then at most two terms in Pkf are nonzero for k large
negative and each coefficient is < C2k/p''. Hence ||-Pfc/||i,p(R) —+ 0 provided
p1 < oo, i.e. p > 1. This shows that
i,k- (2-2.2)
fcez fc€z jez
in the sense of LP(R) convergence. We see that the Haar representation
holds for functions in LP(R) provided p > 1.
But what happens when p = 1? Well, as is typical for orthogonal decom-
positions, the expansion (2.2.2) cannot be valid. Indeed, each of the func-
tions appearing on the right in (2.2.2) has mean value zero. If g € Ia(K)
has mean value zero and / is an arbitrary function from L\(B.), then

f \ f - g \ > \ [ f - I9 =\f f
JM \JR JR \JR
This means that the sum in (2.2.2) cannot possibly converge in Li(M) to /
unless / has mean value zero.
This phenomenon is typical of decompositions for orthogonal wavelets
ip: They cannot represent all functions in Li(R). However, if ip is smooth
enough, the representation (2.2.2) will hold for the Hardy space Hi(R) used
in place of Za(R), and in fact this representation will then hold for functions
in the Hardy spaces Hp(R) for a certain range of 0 < p < 1 that depends on
the smoothness of ip. We shall not discuss further the behaviour of orthogo-
nal wavelets in Hp spaces but the interested reader can consult Prazier and
Jawerth (1990) for a corresponding theory in a slightly different setting.

2.3. Smoothness spaces

We noted earlier the important fact that wavelet decompositions provide a
description of smoothness spaces in terms of the wavelet coefficients. We
wish to illustrate this point with the Haar wavelets and the Lipschitz spaces
in L P (R), 1 <p < oo.
 WAVELETS 9

The Lipschitz spaces Lip(a,L p (R)) of L P (R), 0 < a < 1, I < p < oo,
consist of all functions / € LP(R) for which
\\f-f(-+h)\\Lpm=O(ha), h^Q.
A semi-norm for this space is provided by

0<h<oc
The relationship between the smoothness of / and the size of its Haar coef-
ficients rests on three fundamental inequalities. The first of these says that,
for a fixed k e Z, the Haar functions Hj,k, j € Z, are Lp(R)-stable. Because
of the disjoint support of the Hjtk, j € Z, stability takes the following partic-
ularly simple form: for any sequence (c(j)) £ £P(Z) and S — £ j e z c(j)Hjtk,
we have
\ . (2.3.1)
Sez
This follows by integrating the identity

jez
The other two inequalities are related to the approximation properties of
Sk and the projectors P*:
(J) 11/ - PkfhpW < 2 • 2 - * a | / | L i p ( a , M R ) ) , 0 < a < 1, 1 < p < oo.
(B) | 5 | L i p ( 1 / p , M R ) ) < 2 - 2 f c / " | | 5 | | M R ) , SeSk(X,Lp(R)), 1 < P < oo.
(2.3.2)
The first of these, often called a Jackson inequality (after similar inequali-
ties established by D. Jackson for polynomial approximation), tells how well
functions from Lip(a,L p (R)) can be approximated by the elements of Sk.
The second inequality is known as a Bernstein inequality because of its sim-
ilarity with the classical Bernstein inequalities for polynomials, established
by S. Bernstein.
We shall prove (2.3.2J) and (2.3.2B) for 1 < p < oo. If / G Dk and
h := \I\ = 2~k, then, for all x € I, we obtain

- (]7[ S-h l/(x) " /(x + s ) | P d s ) ' V•

10 R. A. DEVORE AND B. J. LUCIER

If we raise these last inequalities to the power p, integrate over /, and then
sum over all / € Dk, we obtain

11/ "
which implies the Jackson inequality.
The Jackson inequality can also be proved from more general principles.
Since the Pk have norm 1 on LP(R) and are projectors, we have
11/ - Pkf\\Lp(R) < (1 + limi)dist(/,<S fc ) MR ) < 2dist(/,S*) Lp(R) . (2.3.3)
Thus, the Jackson inequality follows from the fact that functions in
Lip(a,L p (E))
can be approximated by the elements of Sk with an error not exceeding the
right-hand side of (2.3.2J).
To prove the Bernstein inequality, we note that any S = J2jez cU)Xij k m
Sk(X,Lp(R)) has norm:
£ 4fcl = E k0T2-*. (2.3.4)

We fix an h > 0. If h > 2~k (i.e. h'1^ < 2k/P), then

/i- 1 /P||5(.+/ l )-5|U p ( R ) <2^(||5(-+/i)|| M
If h< 2~k then
x €
\S(x + K\- S(x)\ - i °'
Therefore
1/P

+ 2"fe) - 5|| L p ( R ) < 2

With the Jackson and Bernstein inequalities in hand, it is now easy to
show that

l/lLiP(Q)Lp(R)) ~ sup2 fc (- +1 /2-i/p)( £ \{f,HjJt)P) , 0 < a < 1/p.

(2.3.5)
(It will be convenient to use the notation A w B to mean that the two
ratios A/2? and 2?/.<4 of the functions A and B are bounded from above
 WAVELETS 11

independently of the variables; in (2.3.5), independently of / . ) First, from

the Jackson inequality,

If we write Pkf - Pk-if = Ej€z(f>HJ,k-i) Hilk-i and replace \\Pkf -

Pk-if\\Lp(R) by the sum in (2.3.1) (with c(j) = (f,Hjtk-i)), we obtain that
the right-hand side of (2.3.5) does not exceed a multiple of the left.
To reverse this inequality, we fix a value of h and choose n G Z so that
2" n < h < 2 " n + 1 . We write / = £ f e 6 Z wk with wk := (Pk+if - Pkf) and
estimate

+h)- wk\\Lp(&) + ^2 \\wk(• + h) - wk\\Lp(R)•

k>n k<n
(2.3.6)
The first sum does not exceed
(mp2ka\\wk\\LAu)) ( |
\k>n / V k>n
Similarly, using the Bernstein inequality, the second sum does not exceed

k<n k<n

SI
>fe<n »)(

fc<n
H
If we write wk = Ejez(/> j,k) Hj,k and use (2.3.1) to replace ||wfc||£p(R) by

in each of these expressions, and then use the resulting expression in (2.3.6),
we obtain

ll/(- +h)-f\\LpW <C/

which shows that the left-hand side of (2.3.5) does not exceed a multiple of
the right.
The restriction a < 1/p arises because the Haar function is not smooth;
for smoother wavelets, the range of a can be increased.
12 R. A. DEVORE AND B. J. LUCIER

2.4- The fast Haar transform

In numerical applications of the Haar decomposition, one must work with
only a finite number of the functions Hj:k- The choice of which functions
to use is often made as follows. Given a function / € Z/2(R), we choose
a large value of n, compatible with the accuracy we wish to achieve, and
we replace / by Pnf with Pn, as before, the Z,2(R) projector onto <Sn, the
space of piecewise-constant functions in £<2(R) with breakpoints at the dyadic
integers j2~", j € Z. If / has compact support then Pnf is a finite linear
combination of the characteristic functions Xi, I G Dn. If / does not have
compact support, it is necessary to truncate this sum (which is justified
because /R\(_OtO] |/| 2 d i -<• 0, a -> oo).
We can now write
n-l
Pnf = (Pnf ~ Pn-lf) + " • + (Plf - POf) + ft/ = ft/ + £ <?*/• (2"41)
fc=0
which is a finite Haar decomposition. We have started this decomposition
with Pof but we could have equally well started at any other dyadic level.
The fast Haar transform gives an efficient method for finding the coeffi-
cients in the expansions
Qkf = £ d(j, k)HiM, d(j, k) := (f, HjM), (2.4.2)
jez
and
Pkf = £ c(j, k)</>iik, c(j, k) := (f, <(>jyk). (2.4.3)
jez
These coefficients are related to the integrals of / over the intervals Ij k '•=
\j2-k,(j + l)2-k):

c(j,k) = 2k'2 I fdx,

d(j,k) = 2k'2[[ fdx- f fdx).

\Jl2j,k+l Jhj + l,k+l )

Therefore, if the coefficients c(j, k + 1), j e Z, are known, then

fy(2j, k + l) + c(2j + l,k + 1)),

d(j,k) =±(c(2j,k + l)-c(2j + l,k + l)). (2.4.4)

In other words, starting with the known values of c(j,n) at level n, we

can iteratively compute all values d(j, k) and c(J, k) needed for (2.4.1) from
(2.4.4). The computation of the c(j,k) at dyadic levels A; ^ 0 is necessary
 WAVELETS 13

for the recurrence even though we are, in the end, not interested in their
values.
There is a similar formula for reconstructing a function from its Haar
coefficients. Now, suppose that we know the coefficients appearing in (2.4.1),
i.e. the values c(j, 0), j G Z, and d(j, k), j € Z, k = 1 , . . . , n, and we wish
to find c(j,n), i.e. to reconstruct / . For this we need only use the recursive
formulae
c(2j,k + l) =fyc(j,k) + d(j,k)),
c(2j + 1, k + 1) = -L(cO\ A:) - d(j, k)).
More information on the structure of the fast Haar transform can be found
in Section 5.

2.5. Multi-variate Haar functions

There is a simple method to construct multi-variate wavelets from a given
uni-variate wavelet, which, for the Haar wavelets, takes the following form.
Let 4>o '•= <t> — X[o,i] a n d <(>i '•= i> = H and let V denote the set of vertices
of the cube O := [0, l]d. For each v = ( « i , . . . , Vd) in V and x = (xi,..., Xd)
from Rd, we let tpv(x) :— nf=i0«i( x i)- The functions tpv are piecewise
constant, taking the values ±1 on the d-tants of fi. The set
* := {t/>v | v 6 V, v^O}
is the set of multi-dimensional Haar functions; there are 2d—1 of them. They
generate by dilation and translation an orthonormal basis for L2(M.d). That
is, the collection of functions 2M/2ipv(2k- - j), j € Zd, k 6 Z, v € V \ {0},
forms a complete orthonormal basis for Z/2(Rd).
Another way to view the multi-dimensional Haar functions is to consider
the shift-invariant space S of piecewise-constant functions on the dyadic
cubes of unit length in Rd. A basis for S is provided by the shifts of X[o,i]d.
Note that the space S is the tensor product of the uni-variate spaces of
piecewise-constant functions with integer breakpoints. The collection of all
shifts of the Haar functions tpv £ \& forms an orthonormal basis for the space
W := S1 Q S°. Properties of the multi-variate Haar wavelets follow from
the uni-variate Haar function. For example, there is a fast Haar transform
and a characterization of smoothness spaces in terms of Haar coefficients.
We leave the formulation of these properties to the reader.

3. The construction of wavelets

3.1. Overview
We turn now to the construction of smoother orthogonal wavelets. Almost
all constructions of orthogonal wavelets begin by using multi-resolution,
14 R . A . DEVORE AND B . J. LUCIER

which was introduced by Mallat (1989) (an interesting exception, presented

by Stromberg (1981), apparently gave the first smooth orthogonal wavelets).
W e begin with a brief overview of multi-resolution that we will expand on
in later sections.
d
Let cp 6 L 2 ( M ) and let <S : = S(</>) be the shift-invariant subspace of
d d
L2(M. ) generated by </>. That is, S(4>) is the L 2 ( R ) closure of finite linear
d
combinations of <f> and its shifts <f>( • + j), j E Z . By dilation, we form the
scale of spaces
k k
S : = {S(2 • ) | S € <S}. (3.1.1)
k k d
Then S is invariant under dyadic shifts j2~ , j € Z . In the construction
of Haar functions, we had d = 1, and «S was the space of piecewise-constant
functions with integer breakpoints. That is, <S = S(<f>) with <j> : = X : = X[o,i]-
Other examples for the reader to keep in mind, which result in smoother
wavelets, are to take for S the space of cardinal spline functions of order r
in L2 ( R ) . A cardinal spline is a piecewise polynomial function defined on
R, of local degree < r, that has breakpoints at the integers and has global
r 2
smoothness C ~ . Then <S = S(N ) with i V the (nonzero) cardinal B-spline r r

that has knots at 0 , 1 , . . . , r. These B-splines are easiest to define recursively:

Ni := X and N : = i V - i * Ni, with the usual operation of convolution
r r

f*g(x) := / f(x-y)g(y)dy.

1
For example, is a hat function, ./V3 a C piecewise quadratic, and so on. In
the multi-variate case, the primary examples to keep in mind are the tensor
product of uni-variate B-splines: N(x) : = N(x\,... ,Xd) •= N(xi) • • • N(xd),
and the box splines, which will be introduced and discussed later.
k
Multi-resolution begins with certain assumptions on the scale of spaces <S
and shows under these assumptions how to construct an orthogonal wavelet
if) from the generating function (f>. The usual assumptions are:
fc k +
(i) «S C S \ k€Z;
d
(ii) [ J ^ = L (R ); 2

5
(hi) D * = {°};
d
(iv) {cf)( - — j)}j d eZ forms an L2(M )-stable basis for <S. (3.1.2)

W e have already seen the role of (ii) and (iii) in the context of Haar decom­
positions. The assumption (iv) means that there exist positive constants C\
and C2 such that each S € S has a unique representation

(i) S = £ c0>(- - j ) , and

d
jez
 WAVELETS 15

2
(ii) Ci||5|| L 2 ( R d ) < ( £ |c(j)| ) ' <C ||5||
2 i 2 ( R e l ) . (3.1.3)

d
If </> has L2(M )-stable shifts then it follows by a change of variables that
kd 2 k d k d
for each k € Z, the function 2 / cp(2 • ) has L ( E ) - s t a b l e 2~ Z shifts. W e 2

shall mention later how the assumption ((3.1.2)(iv)) can be weakened.

Assumption ((3.1.2)(i)) is a severe restriction on the underlying func­
fc
tion <p. Because each space <S is obtained from <S by dilation, we see that
1
((3.1.2)(i)) is satisfied if and only if S C S , or, equivalently, if cj> is in
l d d
the space S . From the Li(R )-stability of the set {cf>{ • — j)}j ' this is ez

equivalent to
<t>(x)= $ > 0 X 2 x - j ) (3.1.4)

d
for some sequence (a(j)) 6 £2{Z ). Equation (3.1.4) is called the refinement
equation for 0, since it says that 4> can be expressed as a linear combination
of the scaled functions 0(2- — j), which are at the finer dyadic level. W e
shall discuss the refinement equation in more detail later and for now only
point out that this equation is well known for the B-spline of order r, for
which it takes the form

r+1
N (x)
r = 2- J2( \ )N (2x-j).
r (3.1.5)
v J J
j=o
Because of ((3.1.2)(i)), the wavelet space
1
w — s es°
1 k
is a subspace of S . By dilation, we obtain the scaled wavelet spaces W,
k fc
k € Z. Then, W is orthogonal to <S and

fc
Since C <S for j < k, it follows that Wj and Wk are orthogonal. Prom
this and ((3.1.2)(ii)) and ((3.1.2)(iii)), we obtain
fc
L (R )
2
d
= 0VT . (3.1.7)
fcez
We find wavelets by showing that W is shift invariant and finding its
generators. For example, when d = 1, W is a principal shift-invariant space,
that is it can be generated by one element ip, i.e. W = S(ip). Of course,
there are many such generators tp for W. In the multi-variate case, the space
d
W will be generated by 2 — 1 such functions.
We find an orthogonal wavelet in one dimension by determining a ip whose
shifts form an orthonormal basis for W. Indeed, once such a function ip is
k 2 k
found, the scaled functions iftj^ : = 2 / tp(2 -— j) will then form an orthonor­
mal basis for L,2(M).
16 R. A. DEVORE AND B. J. LUCIER

Generators xp for W whose shifts are not orthogonal are nonorthogonal

wavelets. For example, if ip has shifts that are I^W-stable (but not or-
thonormal), the functions ipj^ := 2k/2ip(2k- — j) form an Zi2(]R)-stable basis
for £<2(]R). While they do not form an orthonormal system, they still possess
orthogonality between levels,

L >i,k$j',k' dx = 0, k ^ k',

which is enough for most applications. After Battle (1987), we call such
functions if; pre-wavelets.
The construction of (uni-variate) orthogonal wavelets introduced by Mal-
lat (1989), begins with a function (j> that has orthonormal shifts (rather than
just L2(R)-stability). Mallat shows that the function
- j), (3.1.8)
jez
with a(j) the refinement coefficients of (3.1.4), is an orthogonal wavelet.
(It is easy to check that ij) is orthogonal to the shifts of (j> by using the
refinement equation (3.1.4).) A construction similar to that of Mallat was
used by Chui and Wang (1990) and Micchelli (1991) to produce pre-wavelets.
In the construction of pre-wavelets, they begin with a function <j> that has
L2(R)-stable shifts (but not necessarily orthonormal shifts). Then a formula
similar to (3.1.8) gives a pre-wavelet ip (see (3.4.15)).
To find generators for the wavelet space W, we shall follow the construc-
tion of de Boor et al. (1991b), which is somewhat different from that of Mal-
lat. We simply take suitable functions r\ in the space S1 and consider their
orthogonal projections Pi] onto the space 5. The error function w := r\ — Pf]
is then an element of W. By choosing appropriate functions r), we obtain
a set of generators for W. In one dimension, only one function is needed
to generate W and any reasonable choice for r) results in such a generator.
The most obvious choices, r) := 0(2 •) or t] := <f>(2 • — 1), lead to the wavelet
(3.1.8) or its pre-wavelet analogue.
If we begin with the orthonormalized shifts of the B-spline <j> = Nr as
the basis for S, the construction of Mallat gives the spline wavelets ip of
Battle-Lemarie (see, e.g., Battle (1987)), which have smoothness Cr~2. The
support of V" is all of R, although ip does decay exponentially at infinity. More
details are given in Section 3.4. If we do not orthonormalize the shifts, we
obtain the spline pre-wavelets of Chui and Wang (1991), which have compact
support (in fact minimal support among all functions in W).
It is a more substantial problem to construct smooth orthogonal wavelets
of compact support and this was an outstanding achievement of Daubechies
(1988) (see Section 3.5). Daubechies' construction depends on finding a
compactly supported function (j> eCr that satisfies the assumptions of multi-
 WAVELETS 17

resolution and has orthonormal shifts. In this way, she was able to apply
Mallat's construction to obtain a compactly supported orthogonal wavelet
V> in Cr.
The construction of multi-variate wavelets by multi-resolution is based on
similar ideas. We want now to find a set of generators * = {tp} for the
wavelet space W. There are typically 2d — 1 functions in \P. This is an
orthogonal wavelet set if the totality of functions ipjtk, j € Zd, fc 6 Z, ^ 6 $ ,
forms an orthonormal basis for L2(Rd). For this to hold, it is sufficient to
have orthogonality between ip( • — j) and ^ ( - — j ' ) , (j,ip) -^ (j\^)- If the
shifts of the functions ip e # form an Z/2(Rd)-stable basis for W, we say this is
a pre-wavelet set. In this case, we shall still have the orthogonality between
levels: VjVfe J- *l>j',k' if fc ^ k'. Sometimes we also require orthogonality
between ip e \P and all of the tp(- — j), j e Zd, tp ^ t/S. Because the
construction of multi-variate wavelets is significantly more complicated and
more poorly understood than the construction of wavelets of one variable,
we shall postpone the discussion of multi-variate wavelets until Section 3.6.
In the following sections, we shall show how to construct wavelets and
pre-wavelets in the setting of multi-resolution. These constructions depend
on a good description of the space S := S(<j>) in terms of Fourier transforms,
which is the topic of the next section.

3.2. Shift-invariant spaces

Because multi-resolution is based on a family of shift-invariant spaces, it
is useful to have in mind the structure of these spaces before proceeding
with the construction of wavelets and pre-wavelets. The structure of shift-
invariant spaces and their application to approximation and wavelet con-
struction were developed in a series of papers by de Boor et al. (1991a,b,c);
much of the material in our presentation is taken from these references.
We recall that a closed subspace S of Zi2(Rd) is shift invariant if S( • +j),
j G Zd, is in S whenever S £ S. We have already encountered the space
S(<p), which is the I/2(Rrf)-closure of finite linear combinations of the shifts of
<f>. We say that such a space is a principal shift-invariant space (in analogy
with principal ideals). More generally, if $ is a finite set of L^i^) functions,
then the space <S($) is the L2(Rd)-closure of finite linear combinations of the
shifts of the functions <f> G $. Of course, a general shift-invariant subspace
of L2(Rrf) need not be finitely generated.
We are interested in describing the space «S($) in terms of its Fourier
transforms. We let

denote the Fourier transform of an L\ (Rrf) function / . The Fourier transform

has a natural extension from I-i(R d )nL2(R d ) to Z/2(Rd) and, more generally,
18 R. A. DEVORE AND B. J. LUCIER

to tempered distributions. We assume that the reader is familiar with the

rudiments of Fourier transform theory.
The Fourier transform of / ( • +1), t G Rd, is etf; we shall use the abbre-
viated notation
et(x):=eixt
for the exponential functions. Now, suppose that the shifts of <j> form an
Z,2(Md)-stable basis for £(</>)• Then from (3.1.3), each S £ S(4>) can be
written as S = E i € Z d c(j)4>( • - j) with (c(j)) G £2(Zd). Therefore,

s(y) = £ cO>-i(v)*(y) = r(y)ky), <y) •.= £ c{j)e^{y). (3.2.1)

j€Z d j€Zd

Here r is an Z,2(Td) function (i.e. of period 2TT in each of the variables

2/1) • • • > yd)- The L2(Rd)-stability of the shifts of <f> can easily be seen to be
equivalent to the statement

IMIL,(T<) « I | S | | L 2 ( R * ) . (3.2.2)
The characterization (3.2.1) allows one to readily decide when a func-
tion is in <S(<£). Even when the shifts of <f> are not I/2(Rd)-stable, one can
characterize S((f>) by (see de Boor et al. (1991a))

«S(0) = {§ = T<j> € L2(Rd) | T is 27r-periodic}. (3.2.3)

By dilation, (3.2.3) gives a characterization of the scaled spaces Sk, S =
S(<l>). For example, the functions in S1 are characterized by 5 = rfj e
L2(TSLd), rf := <f>(2 •), with r a 47r-periodic function.
A similar characterization holds for a finite set $ of generators for a shift-
invariant space <S($). We say that this set provides L2(Kd)-stable shifts if
the totality of all functions <t>( • — j), j € Z d , <j> G $, forms an I<2(Rd)-stable
basis for «S($). In this case, a function S G S($) is described by its Fourier
transform
= E T<
where the functions r^, <j> G $ , are in LilT*) and

It is clear that the values at points congruent modulo 2TT of the Fourier
transform of a function 5 in S(<f>) are related. If we know 4>(x) and S(x),
then, because r has period 2TT, all other values S(x + a), a G 27rZd, are
determined. It is natural to try to remove this redundancy. The following
 WAVELETS 19

bracket product is useful for this purpose. If / and g are in L2(Kd), we define
[/,<?]:= £ f(-+P)9F+0)- (3-2.4)

Then [f,g] is a function in

One particular use of the bracket product is to relate inner products on
Rd to inner products on T d . For example, if f,g € I<2(Kd) and j G Zd, then

(2w)d f f(x)gjx~^])dx = f ejMfWgJyidy = f ej(9)[f,g}(9) <W.

JR<* jRd JT*
(3.2.5)
Thus, these inner products are the Fourier coefficients of [/, g]. In particular,
a function / is orthogonal to all the shifts of g if and only if [f,g] = 0 a.e.,
in which case one obtains that all shifts of / are orthogonal to all the shifts
of g (which also follows directly by a simple change of variables).
Another application of the bracket product is to relate integrals over Rd
to integrals over T d . For example, if 5 = T(j> with T of period 2TT, then
(2«)d\\S\\L2m = ||r[<U] 1/2 || i2 (T')- (3-2.6)
Returning to .^(R'O-stability for a moment, it follows from (3.2.6) and
(3.2.2) that the shifts of <j> are L2(Rd)-stable if and only if C\ < [<j>, 4>] < C2,
a.e., for constants C\,C2 > 0. Also, the shifts of <f> are orthonormal if and
only if [cj>, $\ = 1 a.e. For example, if we begin with a function <j> with
^ shifts, then the function 0* with Fourier transform

(3-2.7)

has orthonormal shifts (this is the standard way to orthogonalize the shifts
of 0). Incidentally, this orthogonalization procedure applies whenever [<£, 0]
vanishes only on a set of measure zero in T d , in particular for any compactly
supported <f>. That is, it is not necessary to assume that 0 has Z«2(IRd)-stable
shifts in order to orthonormalize its shifts
The bracket product is useful in describing projections onto shift-invari-
ant spaces. Let 0 be an arbitrary Li(BLd) function and let P := P<j, denote
the L2(Rd) projector onto the space S(<f>). Then for each / G L2(Rd), Pf is
the best L2(Rd) approximation to / from S(<f>). It was shown in de Boor et
al. (1991a) that
(328)
^=M^- --
Here and later, we use the convention that 0/0 = 0. We note some properties
of (3.2.8). First, [f,<f>] is 27r-periodic and therefore the form of Pf matches
that required by (3.2.3). If <j> has orthonormal shifts, then [4>,$] = 1 a.e.,
20 R. A. DEVORE AND B. J. LUCIER

and in view of (3.2.5), the formula (3.2.8) is the usual one for the L,2(M.d)
projector. If 4>/[4>, (ft] is the Fourier transform of an L,2(Rd) function 7 (this
holds, for example, if <f> has L2(Md)-stable shifts), then

jmd (3.2.9)

as follows from (3.2.5). Whenever <\> has compact support and L2(Md)-stable
shifts, the function 7 decays exponentially.
The bracket product is also useful in the description of the properties
of the shift-invariant spaces «S($) that are generated by a finite set $ of
functions from I<2(Rd). The properties of the generating set $ are contained
in its Gramian
() (3.2.10)

This is a matrix of 27r-periodic functions from Li(T d ). For example, the

shifts of the functions in $ form an orthonormal basis for <S($) if and only
if G($) is the identity matrix a.e. on Td. The generating set $ provides an
£,2(Rd)-stable basis for «S($) if and only if G($) and G ( $ ) - 1 exist and are
a.e. bounded on 1"* with respect to some (and then every) matrix norm. For
proofs, see de Boor et al. (1991c).

3.3. The conditions of multi-resolution

The question arises as to when the conditions (3.1.2) of multi-resolution
are satisfied for a function <j> € L2(M.d). We mention, without proof, two
sufficient conditions on (ft for ((3.1.2)(ii)) and (iii) to hold. Jia and Mic-
chelli (1991) have shown that if the shifts of (ft are I^R^-stable, if 4> satis-
fies the refinement equation (3.1.4) with coefficients (a(j)) in £i(Xd), and if
£iez<* \<t>(x+j)\ is in L 2 (T d ), then ((3.1.2)(ii)) and ((3.1.2)(iii)) are satisfied.
On the other hand, in de Boor et al. (1991b) it is shown that ((3.1.2)(ii))
and ((3.1.2)(iii)) are satisfied whenever <ft G I«2(Rci) satisfies the refinement
condition ((3.1.2)(i)) and, in addition, supp[0, $] = T d ; by this we mean that
[<^, (ft] vanishes only on a set of measure zero. In particular, these conditions
are satisfied whenever <j> has compact support and satisfies the refinement
condition. Since these conditions are satisfied for all functions <f> that we shall
encounter (in fact for all functions <f> that have been considered in wavelet
construction by multi-resolution), it is not necessary to verify separately
((3.1.2)(ii)) and ((3.1.2)(iii))—they automatically hold. We also note that
in the construction of wavelets and pre-wavelets in de Boor et al. (1991b) it
is not necessary to assume that <j> has Z/2(Rd)-stable shifts.
We now discuss the refinement condition ((3.1.2)(i)). In view of the char-
 WAVELETS 21

acterization (3.2.3) of shift-invariant spaces, this condition is equivalent to

$ = At}, V:=<tf2-) (3-3.1)

d
for some 47r-periodic function A. If 4> has L2(K )-stable shifts then this con-
dition becomes the refinement equation (3.1.4) and A(y) = 2jez rf aU)e-j/2
in the sense of L2(2Td) convergence.
It was shown in de Boor et al. (1991b) that one can construct generators
for the wavelet space W even when ((3.1.2)(iv)) does not hold. For example,
this condition can be replaced by the assumption that supp 0 = Rd. We also
note that if [0,0] is nonzero a.e., then we can always find a generator 0* for
S with orthonormal shifts, so the condition ((3.1.2)(ii)) is satisfied for this
generator (and the other conditions of multi-resolution remain the same).
However, the generator 0 , does not satisfy the same refinement equation as
<j> (for example, the refinement equation for </>* may be an infinite sum even if
the equation for <j> is a finite sum) and <j>+ may not have compact support even
if <j> has compact support, so the construction that gives 0* is not completely
satisfactory. Furthermore, we would like to describe the wavelets and pre-
wavelets directly in terms of the original <j>. This is especially the case when
<j> does not have Z/2(Kd)-stable shifts, since then we can say nothing about
the decay of (f>* even when </> has compact support. In the remainder of this
presentation, we shall assume that <j> has I/2(Rd)-stable shifts.

3.4- Constructions of uni-variate wavelets

In this section we restrict our attention to wavelets in one variable, because
multi-resolution is simpler and better understood for a single variable than
for several variables. We suppose that <j> satisfies the assumptions (3.1.2) of
multi-resolution and follow the ideas presented in de Boor et al. (1991b).
Fundamentally, the approach of de Boor et al. (1991b) is quite simple.
We take a function 77 G S1 and consider its error w := r] — P77 of best L,2(M.)
approximation by the elements of S° = S = S(<f>). Here P is the £<2(K)
projector onto <S(<^>) given by (3.2.8). Clearly w € W and we shall show
that with any reasonable choice for T), the function w is & generator of W,
i.e. W = S(w). Thus, because of the characterization (3.2.3) of principal
shift-invariant spaces, we can obtain all other generators for w by operations
on the Fourier transform side. Here are the details.
We take TJ := <j>{2 •), which is clearly in S1. Then, w := rj — Prj is in W
and by virtue of (3.2.8) has Fourier transform

g4U. (3.4.1)
It is convenient to introduce (for a function / € L2(Kd)) the abbreviated
22 R. A. DEVORE AND B. J. LUCIER

notation
f ••=[/, f}1/2, (3.4.2)
since this expression occurs often in wavelet constructions. Another descrip-
tion of / is

/ = ( E i/(-+i)i2)172-
V
ie2,rZ<» '

We see that / is a 27r-periodic function, and if / has compact support then

P is a trigonometric polynomial, because of (3.2.5). The analogue of this
function for half-shifts is
1 2
/ 2N /
/:=( E l/(-+i)| ) . (3-4-3)
V
j£4irZ<< '

which is now a 47r-periodic function. In particular / has orthonormal half-

shifts if and only if / = 2~ d / 2 a.e., and L2(Rd)-stable half-shifts if and only
if C\ < f < C2, a.e., for constants C\, C2 > 0.
We return now to the construction of wavelets. We can multiply w by
any 27r-periodic function, and as long as the resulting function is in £2(1*),
it will be the Fourier transform of a function in W. We multiply (3.4.1) by
<^2, which clears the denominator. The result is the function WQ with Fourier
transform
j2 (3.4.4)
We note that WQ has compact support whenever <f> does.
We can calculate the bracket products appearing in (3.4.4) by using the
refinement relation <f> = Afj (see (3.3.1)) with A a 47r-periodic function. For
example, to calculate 0,
~4>2 =

= E j64irZ

Similarly, [fj,4>] = Al? + A( • + 2TT)^ 2 ( • + 2 T T ) . Therefore,

= {\A\2fj2 + \A(-+ 27r)| 2 f( • + 2TT) - ~AAfj2

2
^ - +27T)}J)

- A} A(- +27r)f( • + 2ir)fj.

We can make one more simplification in the last representation for too-
 WAVELETS 23

The function \e\/2{A — A( • + 2ir)} is 27r-periodic. Therefore, dividing by

this function, we obtain the function
i> := 2e_ 1/2 A( • + 27r)rf( • + 2^)17. (3.4.5)
It is easy to see (and is shown in (3.4.14)) that if) is in I<2(K). It follows,
therefore, that ip is in W and S(ip) C W. The following argument shows
that we really have S(i>) = W.
If we replace 77 by 771 := rj( • — 5) (which is also in S1) and repeat the
previous construction, in place of wo we obtain the function wi whose Fourier
transform is
wi = e_ 1 / 2 {A( • + 2n) + A} A(- + 2n)fj2(• + 2n)fj.
Hence, dividing by A{ • + 2TT) + A (which is 27r-periodic), we arrive at the
same function ip. The importance of this fact is that we can reverse these
two processes. In other words we can multiply ip by a 27r-periodic function
and obtain r) — Pr\ (respectively f\\ — Prji). Hence, both of these functions
are in S(i>). Since Pt] is in S(<f>), TJ = Pi) + (r) - Prj) is in S(<j>) + S(r/>).
Similarly, r/i is in this space. Since the full shifts of r\ and 771 generate «S1(0),
we must have W = S(ip). This confirms our earlier remark that W is a
principal shift-invariant space. Since we can obtain tj} from w and WQ by
multiplying by 27r-periodic functions, both w and u>o are also generators of
W.
We consider some examples that show that ip is the (pre)wavelet con-
structed by various authors.
Orthogonal wavelets To obtain orthogonal wavelets, Mallat (1989) begins
with a function <f> that satisfies the assumptions (3.1.2) of multi-resolution
and whose shifts are orthonormal. This is equivalent to (f> = 1 a.e., and (by
a change of variables) to the half-shifts of \/2»7 being orthonormal, i.e. to
fj = \ a.e. When this is used in (3.4.5), we obtain
4> = e_1/2A( • + 27r)77, (3.4.6)
which is the orthogonal wavelet of Mallat. To see that the shifts of ip are
orthonormal, one simply computes
i>2 = \A(- +2 7 r)| 2 rf + |A| 2 rf (• +2TT) = \{\A( • +27r)| 2 + |^| 2 } = 1, (3.4.7)
where the last equality follows from the identity
I = tf = l2 + i2( .+27r) = |A|2rf+ |A(-+27r)|2rf(-+27r)
2
I + |A|2}. (3.4.8)
The Fourier transform identity (3.4.6) is equivalent to the identity (3.1.8).
We note that from the orthogonal wavelet ip of (3.1.8) (respectively (3.4.6)),
24 R. A. DEVORE AND B. J. LUCIER

we obtain all other orthogonal wavelets in W by multiplying ^ by a 2TT-

periodic function r of unit modulus. Indeed, we know that any element
w € W satisfies w — Tif> with r G £2(1"). To have [w,w] = 1 a.e., the
function r must satisfy |r(y)| = 1 a.e. in T.
As an example, we consider the cardinal B-spline Nr of order r. To ob-
tain orthogonal wavelets by Mallat's construction, one need only manipulate
various Laurent series. First, one orthogonalizes the shifts of Nr- This gives
the spline <j> = Mr whose Fourier transform is

j> = &:=%-. (3.4.9)

It is easy to compute the coefficients in the expansion

#r2 = 5 > 0 > - ; - (3.4.10)

In fact, we know from (3.2.5) that this is a trigonometric polynomial whose

coefficients are

a(j) = / Nr(x - j)Nr{x) dx= [ Nr(r + j - x)Nr(x) dx

= [Nr*Nr](j + r) = N2r(j + r), j G Z,
because Nr is symmetric about its midpoint.
The polynomial P2r(z) '•= zr Hjez a(j)z~^ *s ^^ e Euler-Frobenius polyno-
mial of order 2r, which plays a prominent role in cardinal spline interpolation
(see Schoenberg (1973)). It is well known that /?2r has no zeros on \z\ — 1.
Hence, the reciprocal l/p2r is analytic in a nontrivial annulus that contains
the unit circle in its interior. One can easily find the coefficients of reciprocals
and square roots of Laurent series inductively. By finding the coefficients of
P2r , we obtain the coefficients f3(j) appearing in the expansion

4>(x) = K(x) = Yl P(J)*r(x - j). (3.4.11)

j€Z

Because piT has no zeros on \z\ = 1, we conclude that the coefficients /3(j)
decrease exponentially. The spline Nr together with its shifts form an or-
thonormal basis for the cardinal spline space S(Nr). They are sometimes
referred to as the Franklin basis for S(Nr).
Now that we have the spline <f> := Mr in hand, we can obtain the spline
wavelet ip = Af* of Battle-Lemarie' (Battle, 1987) from formula (3.1.8). For
this, we need to find the refinement equation for (j>. We begin with the
refinement equation (3.1.5) for the B-spline Nr, which we write in terms of
 WAVELETS 25

Fourier transforms as Nr = Aofjo with 770 := Nr(2 •) and

a 4?r-periodic trigonometric polynomial. It follows that

$ = Af), r,:=K{2'), A(y) = - ^ = Nr(y/2)N-\y)A0(y).

(3.4.12)
In terms of the B-spline Nr, this gives
4>(V) = e-1/2(y)A(y + 2n)f)(y) = \e_l/2{y)A{y + 2n)N-1(y/2)Nr(y/2).
(3.4.13)
In other words, to find the orthogonal spline wavelet ip of (3.4.13), we need to
multiply out the various Laurent expansions making up A(- + 2n)N~1( • /2).
This gives the coefficients 7(3), j G Z, in the representation

iez
We emphasize that each of the Laurent series converges in an annulus con-
taining the unit circle. This means that the coefficients 7(3) converge expo-
nentially to zero when j —* ±00.
Pre-wavelets For the construction of pre-wavelets, we do not assume that
the shifts of <f> are orthonormal, but only that they are I«2(M)-stable, i.e. we
assume ((3.1.2)(iv)). Then, it is easy to see that the function tp defined by
(3.4.5) is a pre-wavelet. Indeed, we already know that tp is a generator for
W and it is enough to check that it has Z,2(IR)-stable shifts. For this, we
compute

^ = \A{ • + 27r)|2if (• + 27r)if + \A\2fjA?j2(• + 2ir). (3.4.14)

Since the shifts of <\> are L2(K)-stable, so are the half-shifts of 77. This means
that C\ < f\ < C2 for constants C\,C2 > 0. Moreover, the formula

shows that C\ < \A\2 + \A(- + 2TT)|2 < C2, again for positive constants
C\,C2- Combining this information with (3.4.14) shows that ip is bounded
above and below by positive constants, so that ip has L2(K)-stable shifts.
This also shows that ift is in I ^ W - The pre-wavelet ip was introduced by
Chui and Wang (1991) and independently by Micchelli (1991).
We can also find a direct representation for ip in terms of the shifts of
<j>(2 •). For this we need the Fourier coefficients fj,(j) (of e-j/2) f° r the 4?r-
26 R. A. DEVORE AND B. J. LUCIER

periodic function 2Afj:

[
JR
Using this in (3.4.5), we find that
- i), Mi) := / 0 W ( 2 z + j) dx. (3.4.15)

If <£ has compact support, then clearly ip also has compact support. Chui
and Wang (1990) posed the interesting question as to whether ^ has the
smallest support among all the elements in W, to which they gave the fol-
lowing answer. We assume that A is a polynomial, i.e. that <t> satisfies a
finite refinement equation. Next, we note that because W C S1, any w € W
is represented as
(3.4.16)
with B of period 4n. If B = J ^ i m 6(j)e_,y2 is a Laurent polynomial with
6(m)6(M) ^ 0, then to has compact support, and we define the length of B
to be M — m. We know that there are nonzero polynomials B that satisfy
(3.4.16) for some w because r? is a polynomial (since r} has compact support)
and (3.4.5) implies that for Bo := Afj , w is the pre-wavelet ip G W.
BQ may not have minimal length among all such polynomials; however,
because it may be possible to cancel certain symmetric factors from So- To
see this, we write Bo(y) = eM(y/2)P(e~iy^2) with P an algebraic polynomial,
and we let Q(z2) := FIAC-2 ~ ^)> with the product taken over all A with A
and —A both zeros of P. Then, the factorization P(z) = Q{z2)Pif{z) gives
the factorization Bo(y) = r(y)fi*(y) with r a trigonometric polynomial of
period 2TT that does not vanish. Therefore, the function rp* with Fourier
transform
My) = e_1/2(y)B*(y + 2ir)fj(y), B*(y) := T-1(y)B0(y), (3.4.17)
is in W and has smaller length than BQ. A simple argument (which we
do not give) shows that Bm has smallest length. For most pre-wavelets of
interest, B* = BQ.
The problem of finding a wavelet w in the form (3.4.16) with B a polyno-
mial of minimal length, which is solved by w = tp», is not always equivalent
to finding the wavelet with minimal support; here the word 'support' means
the interval of smallest length outside of which w vanishes identically. In
general, there may be wavelets w of compact support that can be repre-
sented by (3.4.16) with B not a polynomial. However, Ben-Artzi and Ron
(1990) show that this is impossible whenever the following property holds:
 WAVELETS 27

Fig. 2. The Chui-Wang spline pre-wavelet for r = 4, which has support

[—3,4]. The vertical scale is stretched by a factor of eight.

The linear combination Sjgz 7C?')0(' ~J) (which converges pointwise, since
(j> has compact support) is identically zero if and only if all the coefficients
^(j) are 0. Under these assumptions, the wavelet i/>* has minimal support
(see de Boor et al. (1991b) for details).
For a pre-wavelet ip, we have the wavelet decomposition

Cj,k(f) ••= I flTk, (3-4.18)

where 7 has Fourier transform 7 = V>/[^,i/>]. This follows from the repre-
sentation (3.2.9) for the projector P from Z/2(R) onto W. It is useful to note
that when ip has compact support, the function 7 will generally not have
compact support because of the division by the bracket product [ip, •$[. Thus,
there is in some sense a trade-off between the simplicity of the pre-wavelet
and the complexity of the coefficient functional.
We consider the following important example. Let (j> := Nr be the cardinal
B-spline of order r, which is known to have linearly independent shifts.
Then, the function tjj in (3.4.15) is a spline function with compact support.
It is easy to see that Afj has no symmetric zeros so that ip has minimal
support. We note also that it is shown in de Boor et al. (1991b) that
the shifts of tp are themselves linearly independent. From formula (3.4.15),
we see that if> is supported on [1 — r,r]. Up to a shift, the spline ip is the
minimally supported spline pre-wavelet of Chui and Wang (1991); see Figure
2.

3.5. Daubechies' compactly supported wavelets

The orthogonal spline wavelets of Section 3.4, which decay exponentially
at infinity, can be chosen to have any specified finite order of smoothness.
It is natural to ask whether orthogonal wavelets can be constructed that
have both any specified finite order of smoothness and compact support. A
28 R. A. DEVORE AND B. J. LUCIER

celebrated construction of Daubechies (1988) leads to such wavelets, which

are frequently used in numerical applications. Space prohibits us from giving
all the details of Daubechies' construction, but the following discussion will
outline the basic ideas.
To construct a compactly supported wavelet with a prescribed smoothness
order r and compact support, one finds a special finite sequence (a(j)) such
that the refinement equation (3.1.4) has a solution <f> € Cr with orthogonal
shifts. The orthogonal wavelet if) of (3.4.5) will then obviously have com-
pact support and the same smoothness. Before we begin, it is necessary to
understand which properties of the sequence (a(j)) guarantee the existence
of a function <f> with the desired properties, i.e. we need to understand the
nature of solutions to the refinement equation (3.1.4). This has been stud-
ied in another context, namely in subdivision algorithms for computer aided
geometric design (see, for example, the paper by Cavaretta et al. (1991)
for a discussion of subdivision). As was pointed out by Dahmen and Mic-
chelli (1990), it is possible to derive part of Daubechies' construction from
the subdivision approach. However, we shall describe Daubechies' original
construction.
Let r be a nonnegative integer that corresponds to the desired order of
smoothness, and let {a(j)) with a(j) = 0, \j\ > m, and a(m) ^ 0, be the
sequence of the refinement equation (3.1.4) for the function <f> we want to
construct. The sequence (o(j)) and the Fourier transform of <f> are related
by
, m
<f>(y) = A(y/2)$(y/2), A(y) := - £ a(j)e~^. (3.5.1)
j=-m

Here we use a slightly different normalization for the refinement function

(A(y) = ^A(2y)). If 4> is continuous at 0 and $(0) = 1, we can, at least in
a formal sense, write
0(y) = lim Ak(y) (3.5.2)
fe—K»

where

Ak{y):=f[A(y/V). (3.5.3)
3=1

We note that A1(y) := Ak(2ky) is a trigonometric polynomial of degree

(2k — \)m. The key to Daubechies' construction is to impose conditions on
A (which are therefore conditions on the sequence (a(j))) that not only make
(3.5.2) rigorous but also guarantee that the function <j> defined by (3.5.2) has
the desired smoothness and has orthonormal shifts.
We first note that if the shifts of (f> are orthonormal then, as was shown
 WAVELETS 29

in (3.4.8),
l-4(y)|2 + \A(y + TT)|2 = 1, y € T. (3.5.4)
The converse to this is almost true. Namely, Daubechies' construction shows
that (3.5.4) together with some mild assumptions (related to the convergence
in (3.5.3)) imply the orthonormality of the shifts of <j>. For this, the following
identities, which follow from (3.5.4) by induction, are useful:

£ \A*k(y + J2~k2n)\2 = 1, k = 1,2.... (3.5.5)

i=o
We next want to see what properties of A guarantee smoothness for </>.
The starting point is the following observation. If / R <j>{x) da; ^ 0, then
integrating the refinement equation (3.1.4) gives 52j€Za(j) = 2. Hence,
A(0) = 1 and A(K) = 0. We can therefore write
MV) = (1 + e iy ) iV a(y), ||a|| Loo(T) = 2~\ a(0) = 2 " " , (3.5.6)
for a suitable integer N > 0, a real number 9, and a function a.
By carefully estimating the partial products Ak, it can be shown that
whenever A satisfies (3.5.6) for some 9 > | , the product (3.5.2) converges
to a function in L2O&) that decays like \x\~0 as |x| —» oo. The limit function
is the Fourier transform of the solution <f> to the refinement equation (3.5.1).
We see that the larger we can make 9 in (3.5.6), the smoother <j> is. For
example, if 9 > r + 1, then <f> is in Cr.
What is the role of the integer N in (3.5.6)? Practically, one must increase
N to find a function a(y) that satisfies (3.5.6) for large 9. In addition, the
local approximation properties of the spaces Sk{<j>) are determined by N;
see Section 5.
Once it is shown that there is a function (j> that satisfies the refinement
equation for the given sequence a(j), it remains to show that <p has com-
pact support and orthonormal shifts. Here the arguments have the same
character as those used to analyse subdivision algorithms for the graphical
display of curves and surfaces. Assume that A satisfies (3.5.4) and (3.5.6)
for some 9 > \ and let X denote the characteristic function of [— | , \). Then
X{y) = (siny/2)/(j//2). We define <f>k to be the function whose Fourier
transform is <f>k(y) := Ak(y)x{2~ky). It can then be shown that

iy -» 0, fc - K J O . (3.5.7)

If -Alt = SjO*(j, fc)e_j, then 52ja*(j,k)x{x — j) has Fourier transform

4k(»)*(y) and My) = Al(2-ky)X(2-ky). Therefore,
4>k{x) = £ a*(j, k)2kX(2kx - j). (3.5.8)
jez
30 R. A. D E V O R E AND B. J. LUCIER

Since the coefficients a*(j, k) of A*k axe 0 for | j | > (2* — l ) m , we obtain that
4>k is supported in [—m,m]. Letting k —> oo, we obtain from (3.5.7) that <j>
is also supported on [ - m , m ] . From (3.5.8), (3.5.5), and the orthonormality
of the functions 2k/2x{2k- - j), j 6 Z, we have

<t>k{x)<t>k{x-t)dx = 2k

Here the last equality follows by expanding the identity (3.5.5). Letting
A; —> oo, we obtain that {<f>{ • — j)}jez is an orthonormal system.
This outline shows that a Cr, compactly supported function <j> with or-
thonormal shifts exists if (3.5.4) and (3.5.6) hold for a sequence (a(j)) and
two numbers N and 6 > r + 1. The following arguments show that such
sequences exist.
We look for an A of the form (3.5.6) with a a trigonometric polynomial
with real coefficients. Then, |a(y)| 2 = a(y)a(y) = a(y)a(—y) is an even
trigonometric polynomial, and
\a(y)\2 = T(cosy) = T(l - 2sin 2 y/2) = R(sin2 y/2)
with R an algebraic polynomial. The identity (3.5.4) now becomes
(cos22//2)JVJR(sin2j//2) + (sin2 i//2)Arii(cos2 y/2) = 2~2N.
With t := sin2 y/2, we have
(1 - t)NR(t) + tNR(l -t) = 2~2N. (3.5.9)
Therefore, to find A, we must find an algebraic polynomial R that satisfies
(3.5.9). It is easy to see that the degree of R must be at least N — 1. We
can find R of this degree by writing R in the Bernstein form

fc=o
Then, (3.5.9) becomes

fc=O \ / fc=O

k
k=o \ )
We see that
(2N-1\
\, :
r,- 1 )'
= 2o-2JV \ fc / u —n i ^r 1
 WAVELETS 31

satisfies (3.5.10), and we denote the polynomial with these coefficients by

RN-
It is important to observe that !?#(£) is nonnegative for 0 < t < 1, because
we wish to show that there is a trigonometric polynomials a(y) such that
R,N(sin2y/2) = \a(y)\2, i.e. we somehow have to take a 'square root' of
RN- For this, we use the classical theorem of Fejer-Riesz (see, for example,
Karlin and Studden (1966) p. 185) that says that if R is nonnegative on
[0,1], then i?(sin 2 y/2) = |a(y)| 2 for some trigonometric polynomial a with
real coefficients and of the same degree as R. We let CXN be the trigonometric
polynomial corresponding to RN.
We now set AN(V) '•= (1 +eiv)N(>N(y) and note that AN satisfies (3.5.4)
because RN satisfies (3.5.9). Therefore, the function <f> defined via the limit
process (3.5.2) has compact support and orthonormal shifts. The corre-
sponding orthogonal wavelet ip =: V2N defined by (3.1.8) for the refinement
coefficients (a(j)) is an orthogonal wavelet with compact support. It is easy
to show that V2N is supported in [— (N — 1), JV].
The question now is what is the smoothness of V2N- Here the matter
can become somewhat technical (see Daubechies (1988) and Meyer (1987)).
However, the following 'poor man's' argument based on Stirling's formula
at least shows that given any integer r, if we choose N sufficiently large, the
orthogonal wavelet 2?2JV will have smoothness C r .
Because the Bernstein coefficients of RN are monotonic, it follows that
RN is increasing on [0,1]. Therefore, maxo<t<i-RJV(*) = -RAT(I) = XN-I =
2-2N(2N~1). Therefore, Hajvlll^) is bounded by

\ N )

by Stirling's formula. We see that given any value of 9 > 0, we can choose
N large enough so that (3.5.6) is satisfied for that 0, and the function 4>
satisfies \4>(x)\ < C(\ + \x\)~6. Hence, for any r < 9 — 1, (p, and hence VIN,
is in C.
For N = 1, the Daubechies construction gives <j> = X[o,i] and T>2 is the
Haar function. For N = 2, the polynomial R2(t) = (1 + 2tj/16 and

A2(y) = (1 + e*)2 (^±± _ ^ l i e - ^ = (1 + J»)2a2(y).

Then a2 (y) satisfies |«2(l/)| < N/3/4 < 2" 1 . Therefore, the function <f>
and the wavelet Vi := ip corresponding to this choice is continuous. (See
Figure 3 for a graph of <f> and tjj.) A finer argument shows that V4 is in
Lip(.55,Loo(R)). The reader can consult Daubechies (1988) for a table of
32 R. A. D E V O R E AND B. J. LUCIER

Fig. 3. The function <j> and the Daubechies wavelet rjj = V2N when JV = 2.

the refinement coefficients of Z?2JV for other values of N and a more precise
discussion of the smoothness of Z>2N in Loo(R)-

3.6. Multi-variate wavelets

There are two approaches to the construction of multi-variate wavelets for
I»2(Rd). The first, the tensor product approach, we now briefly describe.
In this section, V will denote the set of vertices of the cube [0,1]d and
V := V \ {0}. Let <j> be a uni-variate function satisfying the conditions
(3.1.2) of multi-resolution and let ip be an orthogonal wavelet obtained from
(f>. For <f>o := <f>, <j>\ := rp, the collection ^ of functions
ipv(xi,...,xd) := veV, (3.6.1)
generates, by dilation and translation, a complete orthonormal system for
L,2(M.d). More precisely, the collection of functions ipj,k,v '•= 2k^2ipv(2k-—j),
j € Zd, k € Z, v G V, forms a complete orthonormal system for d
d
each / € I/2(M. ) has the series representation
/ = (/• ^.*. (3.6.2)

in the sense of convergence in L,2(Rd). This construction also applies to

pre-wavelets, thereby yielding a stable basis for Z,2(Rd).
Another view of the tensor product wavelets is the following. We let <S
be the space generated by the shifts of the function x i—• <j>{x{) • • -0(xd).
Then, the wavelets t(>v are generators for the wavelet space W := S1 Q S°.
The second way to construct multi-variate wavelets uses multi-resolution
in several dimensions. We let </> be a function in L2(Kd) that satisfies the
 WAVELETS 33

conditions (3.1.2) of multi-resolution for S := S((j>), and we seek a set \Jr

of generators for the wavelet space W := S1 0 5°. For example, if we
want an orthonormal wavelet basis for L2(Kd), we would seek \& such that
the totality of functions ip( • — j), j e Zd, ip € * , forms an orthonormal
basis for W. By dilation and translation, we would obtain the collection
of functions ipjtk := 2kd/2ip(2k- - j), ip G * , j G Zd, k G Z, which together
form an orthonormal basis for L,2(M.d). Each function / in Li2(M.d) has the
representation
EEE(/>fe)fe- (3-6-3)
fc€Z jgz<*

For a pre-wavelet set $ , we would require L2(Rd)-stability in place of or-

thogonality. Sometimes, we might require additionally that the shifts of ip
and those of ip are orthogonal whenever ip and %j) are different functions in

Constructing orthogonal wavelets and pre-wavelets by this second ap-

proach is complicated by the fact that there does not seem to be a straight-
forward way to choose a canonical orthogonal wavelet set from the many
possible wavelet sets it.
The book of Meyer (1990) contains first results on the construction of
multi-variate wavelet sets by the second approach. This was expanded upon
in the paper of Jia and Micchelli (1991). These treatments are not always
constructive; for example, the latter paper employs in some contexts the
Quillen-Suslin theorem from commutative algebra. Several papers (Riemen-
schneider and Shen, 1991, 1992; Chui et al., 1991; Lorentz and Madych,
1991) treat special cases, such as the construction of orthogonal wavelet
and pre-wavelet sets when 4> is taken to be a box spline. The paper of
Riemenschneider and Shen (1991) is particularly important, since it gives a
constructive approach that applies in two and three dimensions to a wide
class of functions 4>.
We shall follow the approach of de Boor et al. (1991b), which is based on
the structure of shift-invariant spaces. This approach immediately gives a
generating set for W, which can then be exploited to find orthogonal wavelet
and pre-wavelet sets. De Boor et al. start with a function <j> that satisfies
the refinement relation ((3.1.2) (i)) and whose Fourier transform satisfies
supp<^ = Rd. It is not necessary in this approach to assume ((3.1.2)(ii))
and ((3.1.2)(iii)) - they follow automatically. It is also not necessary to
assume ((3.1.2)(iv)). In particular, this approach applies to any compactly
supported <j>. To simplify our discussion, we shall assume in addition to
((3.1.2)(i)) that <j> has compact support and that the shifts of <f> are L2(Kd)-
stable; we refer the reader to de Boor et al. (1991b) for a discussion of the
more general theory.
The usual starting point for the construction of multi-variate wavelets is
34 R. A. DEVORE AND B. J. LUCIER

the fact that the dilated space S1 of S := S(<f>) is generated by the half-
shifts of r) := <t>(2-), and therefore also by the full shifts of the functions
r/v := rf(- — v/2), v € V. The assumption that supp <j> — Rd is important
because it implies that the set $ := {<f>v := (j>(x — v/2) \ v € V} is also a
generating set for S1, i.e. S1 = <S($). $ is more useful than {r)v | V G V}
as a generating set because $ contains a function that is in 5 ° , namely <j>.
In analogy with the uni-variate construction, we see that with P the Z/2(Kd)
projector onto S, the functions <f>v — P<f>v, v € V, form a generating set for
W. Prom (3.2.8), we calculate the Fourier transforms of these functions and
multiply them by [<£, <j>] to obtain the functions wv, with Fourier transform
tiv:=[$,4]+v-[iv,4>]t, veV. (3.6.4)
The set W := {wv | v 6 V'} is another generating set for W. We note
that because we assume (f> has compact support, the two bracket products
appearing in (3.6.4) are trigonometric polynomials, and hence the functions
wv also have compact support.
The set TW, with T = (TVtVi)VtV>^v' a matrix of 27r-periodic functions, is
another generating set for W if det(T) ^ 0 a.e.
It is easy to find an orthogonal wavelet set by this approach. Because the
functions in W have compact support, the Gramian matrix ([wv, uv])u,u'ev
has trigonometric polynomials as its entries. Since this matrix is symmetric
and positive semi-definite, its determinant is nonzero a.e. We can use Gauss
elimination (Cholesky factorization) without division or pivoting to diago-
nalize G(W). That is, we can find a (symmetric) matrix T = (rv<vi)VtVi&v' of
trigonometric polynomials such that W* := TW has Gramian G(TW*) =
TG(W)T* that is a diagonal matrix with trigonometric polynomial entries.
If w* are the functions in VV*, then the functions w^* with Fourier trans-
forms w** := w*/[w*,w*\1/2, v 6 V, have shifts that form an orthonormal
basis for W. Indeed,

which shows that the new set of generators W** has the identity matrix as
its Gramian.
The disadvantage of the orthogonal wavelet set W** is that usually we
can say nothing about the decay of the functions w**, since this construction
may involve division by trigonometric polynomials that have zeros. How-
ever, when (j> has Zr2(lRd)-stable half-shifts, this construction can be modified
to give an orthogonal wavelet set whose elements decay exponentially (see
de Boor et al. (1991b)).
While the assumption that the half-shifts of <f> are L2(Kd)-stable is often
not realistic, we shall assume it a little longer in order to introduce some new
ideas that can later be modified to drop the stability assumption. Under the
 WAVELETS 35

half-shift stability assumption, we have that ^ is a trigonometric polynomial

of period 4TT that has no zeros. Therefore, 0* := <j>/<j> serves to define an
L 2 (E d ) function in 51(^>) that decays exponentially and has orthogonal half-
shifts. Moreover, the function w with Fourier transform w := <j>/4> is also
in <S1((^) and decays exponentially. Therefore, with [•, -}\/2 the bracket
product for half-shifts (which is defined as in (3.2.4) except that the sum is
taken over 47rZd), we have
&«>]i/2 = $*>&•] 1/2 = 1 a.e. (3.6.5)
The Fourier coefficients (with respect to the e_j/ 2 , j G Ir) of [</>, u)]i/2 are
the inner products of <j> with half-shifts of w. Hence, all nontrivial half-
shifts of w are orthogonal to <f>. This means that the functions in Wo :=
{w( • + v/2) | v e V'} are all in W. It is easy to see that they generate W,
that is W = S(Wo).
Thus, in the special case we are considering, W is generated by the non-
trivial half-shifts of a single function w. It is natural to ask whether this
holds in general (i.e. when we do not assume stability of half-shifts). To see
that this is indeed true, we modify the argument in (3.6.5). If we multiply
w by the 27r-periodic function IlAe27rV 4>(' + ^ ) 2 ' the result is a compactly
supported function w* € Z-2(Rd), with Fourier transform

™* := $ I I U • + A ) 2 - (3-6-6)
Ae2irV
We find that
k-+V2- (3-6.7)
Because the right-hand side is 27r-periodic, we deduce that the inner product
of <f> with w*( • — j/2) is zero whenever j = v + 2k with v EV' and k G Z d .
Hence, the functions tw*( • + v/2), v € V, are all in W and it is easy to see
that they are also a generating set for W.
While the nontrivial half-shifts of w* are a generating set for W, they
have the drawback that they usually do not provide an L2(Rd)-stable basis.
The usual approach towards constructing an L2(Kd)-stable basis for W is to
begin with the generating set {r)v | v € V}, t] := <f>{2 •), for 5 X (0). With this
as a starting point, Meyer (1990) III, Section 6 and Jia and Micchelli (1991)
have shown the existence of a set of generators for W consisting of compactly
supported functions whose shifts provide an L2(Kd)-stable basis for W. How-
ever, their proofs are not constructive. In one, two or three dimensions, and
with an additional assumption on the symmetry of <j>, Riemenschneider and
Shen (1991, 1992), have given a constructive proof for the existence of such
a generating set, which we now describe.
36 R. A. DEVORE AND B. J. LUCIER

We begin again with the function w :— (/)/</> . If r = Jljezd cU)e-j/2 i s a

47r-periodic function whose Fourier coefficients c(j) = 0 whenever j e 2Zd,

then the function with Fourier transform TW is in W provided it is in L2(Rd).
The condition on these Fourier coefficients is equivalent to requiring that
£ r( • + v) = 0 a.e. (3.6.8)
v€2nV
In particular, we can use this method to produce functions in W as follows.
We assume that <f> is real valued, has I/2(K<i)-stable shifts, and satisfies the
refinement equation (3.3.1) for a real valued function A. This assumption
on A is a new ingredient; it will be fulfilled for example if </>(—x) = <j>(x).
(More generally, one only needs symmetry about the centre of the support
of 4>.) If a 6 V and va G 2wV, then we claim that the function ipa with
Fourier transform
i>a •= 2 e Q / 2 (Af)2)( • + vQ)fj (3.6.9)
is in W, provided ea/2(va) — —1. Indeed, the refinement equation says that
io = fj/{Ar) ) . We obtain %pa from w by multiplying by Ta := 2ea/2Aij A( • +
va)fj (• + va). The vertices v and v + va (modulo 2TT) contribute values in
(3.6.8) that are negatives of one another. Hence, (3.6.8) is satisfied for r a and
ipa G W. The functions ipa have compact support when A is a polynomial
and <j> has compact support.
We are allowed to make any assignment of a \—> va with ea/2{va) = — 1-
To obtain an L2(Rrf)-stable basis for W, we need a special assignment with
the property that a — /? (modulo 2) is assigned va — vp (modulo 2n). If such
a special assignment can be made, then a simple computation shows that
""'2

with fj, := 2Afj ,

$a,fo]= £ e(a_fi)/2{-+v)n{-+v+va)ii{-+v+V(i)fj2{-+v). (3.6.10)
v€2irV
For example, if the shifts of <f> are orthonormal, then fj = | . In this case,
for a 7^ 0, the terms of the sum in (3.6.10) are negatives of one another for
the two values v and v + (va + v@) (this is where we need to assume that a
special assignment exists), and hence the sum in (3.6.10) is 0. When a = /?,
this sum is
A2{+)~f\+) Yl 4?(-+v) = 4? = l a.e.

Hence, the Gramian of \P := {ipa \ a £ V'} is the identity matrix. We obtain

an orthonormal basis for W in this way.
If we begin with a 0 that has L2(Md)-stable shifts then a slightly more
complicated argument shows that the functions ipa, a G V, are an Z<2(Rd)-
stable basis for W.
 WAVELETS 37

This leaves the question of when we can make an assignment a i—> va of

this special type. Such assignments are possible for d = 1,2,3 but not for
d > 3. For example, when d = 2, we can make the assignments as follows:

(0,0) — 2 * ( 0 , 0 ) , (0,1) — 2 i r ( 0 , l ) ,
(1,0) — 2 i r ( l , l ) , (1,1) — 2 * r ( l , 0 ) .
This construction will give orthogonal wavelet and pre-wavelet sets for
box splines. To illustrate this, we consider briefly the following special box
splines on R2. Let A be the triangulation of R2 obtained from grid lines
x\ = k, X2 = k, and X2 — x\ =fc,k € Z. Let M be the Courant element for
this partition. Thus, M is the piecewise linear function for this triangulation
that takes the value 1 at (0,0) and the value 0 at all other vertices. The
Fourier transform of M is
sin(y 2 /2)\ /sin(( y i + y 2 )/2)\
~^j-){ (W+W)/2 )•
By convolving M with itself, we obtain higher order box splines defined
recursively by M\ :— M and Mr := M * MT-\. Then Mr is a compactly
supported piecewise polynomial (with respect to A) of degree 3r — 2 and
smoothness C 2 r ~ 2 . Since M is real, the box spline Mr satisfies the refinement
identity (3.3.1) with A real. Therefore, if we take <j> := Mr and <S = S(Mr),
the construction of Riemenschneider and Shen applies to give a pre-wavelet
set ty consisting of three compactly supported piecewise polynomials for the
partition A/2. The set * provides an L2(Ed)-stable basis for the wavelet
space W.

4. Fast wavelet transforms

It is easy to compute the coefficients in wavelet decompositions iteratively
with a technique similar to the fast Haar transform. We shall limit our
discussion to Daubechies' orthogonal wavelets with compact support in one
dimension. However, the ideas presented here apply equally well to other or-
thogonal wavelets and to pre-wavelets. We let <p be a real-valued, compactly-
supported function with orthonormal shifts that satisfies the conditions of
multi-resolution and in particular the refinement equation (3.1.4). The func-
tion <j> is real and the refinement coefficients are real and finite in number.
The orthogonal wavelet %j) is then obtained from <f> by (3.1.8).
A numerical application begins with a finite representation of a function
/ as a wavelet sum. This is accomplished by choosing a large value of n,
commensurate with the numerical accuracy desired, and taking an approx-
imation to / of the form
Sn = £ / A n , (4.1)
jGZ
38 R. A. DEVORE AND B. J. LUCIER

with only a finite number of nonzero coefficients fj. The coefficients fj are
obtained from / in some suitable way. For many applications, it suffices to
take fj := f(j2~n). The point j2~n corresponds to the support of $jtn.
Since Sn € Sn, we have
k
Sn = PnSn = PoSn +f f^iPkSn-Pk-iSn) = P0Sn + E £ c(j, 2
k)^{2kX-j).
fc=l fc=0j€Z
(4.2)
We will present an efficient algorithm for computing the coefficients c(j, k)
from the fj and an efficient method to recover Sn from the coefficients c(j, k).
The algorithm presented later has two main features. First, it computes
the coefficients c(j, k) using only fj and the coefficients a(j) of the refinement
equation (3.1.4) for <j>. In other words, one never needs to find a concrete
realization of the functions (j> and ip. Second, the iterative computations are
particularly simple to program and run very quickly—the complexity of the
fast wavelet transform of 2™ coefficients is O(2n); in contrast, the complexity
of the Fast Fourier Transform is O{n2n).
During one step of our algorithm, we wish to find the coefficients of
Pk-iS when S = E? G Z s 0')^i,* *s m <^fc- The coefficients of Pk-iS =
2t€z s 'W0i,fc-i a r e * n e inner product of S with the 4>i,k-i- We therefore
compute, using (3.1.4),

? E «(*)*a+d = -4 E «(i - a)
In other words, the sequence s' := (s'(i)) is obtained from a := (s(i)) by
matrix multiplication:

a' = As, A:={aij), aitj:=-j=a{j~2i), i,j€Z. (4.3)

Let Qk be the projector onto the wavelet space Wk. A similar calculation
tells us how to compute the coefficients t = (t(i)) of the projection Qk-iS =
E i € z *(»)^,fc-i of S € «Sfe onto Wk-V

t = Bs, B := (ftj), fly := ^ 6 ( j - 2i), i,j G Z, (4.4)

where bj := (—l) J a(l — j ) are the coefficients of the wavelet ip given in

(3.1.8).
This gives the following schematic diagram for computing the wavelet
coefficients (c(j,k)):
 WAVELETS 39

A
Sn ~»
\B \B \B (4.5)
Qn-lSn ... QoSn
In other words, to compute the coefficients of Pk-\Sn we apply the matrix
A to the coefficients of PkSn, while to compute those of Qk-iSn we apply
the matrix B to the coefficients of PkSn- The coefficients c(j, k), j € Z, are
the coefficients of QkSn- This is valid because Qk-iPkSn = Qk-iSn and
Pk~lPkSn = Pk-lSn.
Now suppose that we know the coefficients of PoSn and QkSn, k =
0,... ,n — 1. How do we reconstruct 5 n ? We need to rewrite an element
5 <E Sk, S = J^jez sU)^>j,k as an element of Sk+1, S = £ i e z ,
Prom the refinement equation (3.1.4), we find

S = Ys(j)[~Ya{
Lv
j€Z * t€Z
Therefore, we can express the computation of a' from s as multiplication by
the transpose A* of A:

s' = A*s, A* :=«,), alj:=-j=a(i-2j), ijez. (4.6)

A similar calculation tells us how to rewrite a sum 5 = £ j € Z t(i)0i,fc as a

SUm S ^ zZiPZ ^ (v*ri fc-t-1*

t' = B*t, B* :=(/?,*,-),

/?,*,-:= -7=6(1 - 2j), t , j 6 Z . (4.7)
V2
The reconstruction of Sn from QkSn, k = 0 , . . . , n — 1, and Fo^n is then
given schematically by
A* A* A*
PoSn -» P i 5 n - ^ ...-*5n
/ S* / B* (4.8)
• •• Qn-lSn
The matrices A, B, A*, and B* have a small finite number of nonzero
elements in each row, so each of the operations in (4.5) and (4.6) has com-
putational complexity proportional to the number of unknowns.
The reconstruction algorithm can be used to display graphically a finite
wavelet sum S. We choose a large value of n, and use the reconstruction
algorithm to write 5 = ^ € Z s(j, n)^, n < The piecewise linear function with
values s(j, n) at j2~n is an approximation to the graph of S. Such procedures
40 R. A. DEVORE AND B. J. LUCIER

for graphical displays are known as subdivision algorithms in computer aided

geometric design.
The matrices A and B have many remarkable properties summarized by
BA* = 0,
AA* = I, BB* = I,
A* A + B*B = / . (4.9)
The first equation represents the orthogonality between W and «S°, the sec-
ond the fact that the shifts of <f> and ip are orthonormal, and the third the
orthogonal decomposition S1 — S° © W.

5. Smoothness spaces and wavelet coefficients

We have seen in Section 2 that one can determine when a function / €
LP(R) is in a Lipschitz space Lip(a,Lp(M)), 0 < a < 1/p, by examining the
coefficients of the Haar expansion of / . In fact, one can often characterize
membership in general smoothness spaces in terms of the size of coefficients
in general wavelet or pre-wavelet expansions. We do not have the space
to explain in detail how such characterizations are proved, but we shall
outline one approach, based on approximation, that parallels the arguments
in Section 2 about Haar wavelets. A more complete presentation, much along
the lines given here, can be found in the book of Meyer (1990). The article
of Frazier and Jawerth (1990) gives a much more general and extensive
treatment of wavelet-like decompositions from the viewpoint of Littlewood-
Paley theory.
We shall suppose that 4> satisfies the conditions (3.1.2) of multi-resolution.
We also suppose that <f> has compact support. This is not a necessary as-
sumption for the characterizations given below (it can be replaced by suit-
able polynomial decay at infinity) but it will simplify our discussion. We
shall also assume that 1 < p < oo. The arguments that follow can be
modified simply to apply when p — oo; the analysis for p < 1 can also be
developed as shown later, but then it must be carried out in the setting of
the Hardy spaces Hp(Rd).
We fix a value of p and let S := S(<j>, Lp(Rd)) be the Lp(Rd) closure of the
finite linear combination of shifts of (p.
We assume that the shifts of <j> form an Lp(Rd)-stable basis for S. For
functions with compact support, this holds whenever the shifts of <j> form an
Z/2(Md)-stable basis for S(<f>, L2(M.d)) (see Jia and Micchelli (1991)). It follows
that the dilated functions <f>jtk,p := 2k/p</>(2k- - j), j e Zd, form an Lp(Rd)-
k
stable basis of S , for each k € Z. That is, there are constants C\,Ci > 0
such that each S € Sk can be represented as S = J2jezd c0>
 WAVELETS 41

with
\ i/p

n < \\s\\Lp{Rd) < c 2

J
We will now show that the orthogonal projector P from L2(Rd) onto
S(<j>, L,2(Rd)) has a natural extension to a bounded operator from Lp(d
onto 5 . We can represent P as in (3.2.9):

(5-2)

The function 7 € Loo(Rd) decays exponentially at infinity and hence is in

Lq(Rd) for 1 < q < 00. In particular, 7 e Lpi(Rd), and (5.2) serves to define
P on L p (R d ). The compact support of <j> and the exponential decay of 7 then
combine to show that P is bounded on Lp(Rd). By dilation, we find that
the projectors Pf. (which map Lp(Rd) onto Sk) are bounded independently
of it.
The projector Q from L2(Rd) onto the wavelet space W is also represented
in the form (5.2) and has an extension to a bounded operator on Lp(Rd) for
the same reasons as before. We can also derive the boundedness of Q from
the formula Q = P\ — PQ.
Since the Pk are bounded projectors onto Sk, their approximation prop-
erties are determined by the approximation properties of the spaces <Sfc.
Consequently, we want to bound the error of approximation by elements in
Sk of functions in certain smoothness classes. In particular, we are interested
in determining for which spaces Sk it is true that
dist(/,«S fe ) Lp(R , ) < C2-kr\f\wr{Lp{Rd)); (5.3)
here Wr(Lp(Rd)) is the Sobolev space of functions with r (weak) derivatives
in Lp(Rd) with its usual norm and semi-norm. This well-studied problem
originated with the work of Schoenberg (1946), and was later developed by
Strang and Fix (1973) for application to finite elements. Strang and Fix
show that when <j> has compact support, a sufficient condition for (5.3) to
hold is that
0(0) # 0 and £>"0(27ra) = 0, |i/| < r, a G Zd, a ^ O . (5.4)
This condition is also necessary in a certain context; see de Boor et al.
(1991a) and de Boor and Ron (1991) for a history of the Strang-Fix condi-
tions.
Schoenberg (1946) showed that (5.4) guarantees that algebraic polynomi-
als of (total) degree < r are contained locally in the space Sk. This means
that for any compact set Q and any polynomial R with deg(R) < r, there
is an S € S that agrees with R on fi.
42 R. A. DEVORE AND B. J. LUCIER

Fig. 4. The function 0 • <j>{x + 1) + 1 • <j>(x) + 2 • <j>(x - 1), with <l> given by
Daubechies' formula (3.5.1) with N = 2; see Figure 3. Note the
linear segment between x = 1 and a; = 2.

In summary, the approximation properties of Pk are determined by the

largest value of r for which (5.4) is valid. Because we usually know a lot
about the Fourier transform of <j>, the best value of r is easy to determine.
For this r, we have (5.3). For example, for <f> = Nr, the B-spline of order
r, <j>(y) = (1 — e~iy)r/(iy)r satisfies (5.4) for this value of r. Similarly, the
Daubechies wavelets satisfy (5.4) for r = N, with N the integer appearing
in the representation (3.5.6); see Figure 4.
Assume now that there are positive integers r and s such that the following
Jackson and Bernstein inequalities hold:
(j) ||/

(B) |5b s ( L p ( R d ) ) <C-2 f e s ||5|| L p ( R d ) ) SeSk(<t>,Lp(R% l<p<oo.

(5.5)
(Actually, s need not be an integer.) The Jackson inequality is just a refor-
mulation of (5.3), and the largest value of r for which (J) holds is determined
by (5.4). The Bernstein inequality holds if <j> G Ws{Lp(Rd)) and in partic-
ular (since (j> has compact support) whenever <f> is in C. It is enough to
verify (B) for k = 0, since (B) would then follow for general k by rescaling.
The left semi-norm in (B) for S = J2jezd CU)H" ~ J) is bounded by the
£p(Zd) norm of the coefficients (c(j)) j€Z d, which is bounded in turn by the
right-hand side of (B) by using the £-p(E<i)-stability of the shifts of <j>.
Once the Jackson and Bernstein inequalities have been established, we
can invoke a general procedure to characterize smoothness spaces in terms
of wavelet coefficients. To describe these results, we introduce the Besov
 WAVELETS 43

spaces, which are a family of smoothness spaces that depend on three pa-
rameters. We introduce the Besov spaces for only one reason: They are
the spaces that are needed to describe precisely the smoothness of functions
that can be approximated to a given order by wavelets. The following dis-
cussion is meant as a gentle introduction to Besov spaces for the reader who
instinctively dislikes any space that depends on so many parameters.
The Besov space B%(Lp(Rd)), a > 0 and 0 < q,p < oo, is a smooth-
ness subspace of Lp(Rd). The parameter a gives the order of smoothness
in Lp(Rd). The second parameter q gives a finer scaling, which allows us
to make subtle distinctions in smoothness of fixed order a. This second
parameter is necessary in many embedding and approximation theorems.
To define these spaces, we introduce, for h € Rd, the rth difference in the
direction h:

Arh(f,x) := J2(-Vr+J ( \ ) v J
j=o '

Thus, Ah(f,x) := f(x + h) — f(x) is the first difference of / and the other
differences are obtained inductively by a repeated application of A/,. With
these differences, we can define the moduli of smoothness
w r (/, t)p := sup ||AJi(/, • )|| L p ( R d ) , t > 0,

for each r = 1,2,— The rate at which ujr(f,t)p tends to zero gives in-
formation about the smoothness of / in Lp(Rd). For example, the spaces
Lip(a,Lp(Rd)), which we have discussed earlier, are characterized by the
condition wi(/,i) p = O(ta), 0 < a < 1.

The Besov spaces are defined for 0 < a < r and 0 < p, q < oo as the set

of all functions / e Lp(M.d) for which

U J (5 6)
V * (t*\ -
a
Because we allow valuessup* u>(ft)
of p and q oodoes not always
q less than one, this 'norm'
satisfy the triangle inequality, but it is always a quasi-norm: There exists a
constant C such that for all / and g in B%(Lp(Rd)),

Even though the definition of the B^{Lp{Rd)) norm depends on r through

the modulus of smoothness, we have not parametrized the spaces by r, for
two reasons. First, no one can stand spaces that are parametrized by more
than three parameters. Second, it can be shown that all values of r greater
44 R. A. DEVORE AND B. J. LUCIER

than a give rise to equivalent norms, so the set of functions in B™(Lp(Rd))

does not depend on r as long as r > a.
We note that the family of Besov spaces contains both the Lipschitz
spaces Lip(a,L p (R d )) = B^(Lp(Rd)), 0 < a < 1, and the Sobolev spaces
Wa(L2(Rd)) = B%(L2(Rd)), which are frequently denoted by Ha(Rd).
We have mentioned that once Jackson and Bernstein inequalities have
been established, there is a general theory for characterizing membership in
Besov spaces by the decay of wavelet coefficients. This is based on results
from approximation theory described (among other places) in the articles
by DeVore et al. (1991d), DeVore and Popov (1988), and the forthcoming
book by DeVore and Lorentz (1992). Among other things, this theory states
that whenever (5.5) holds, we have that

(5-7)

for 0 < a < min(r, s), 1 < p < oo, and 0 < q < oo.
If * is a wavelet set associated with <j>, then
Qk(f) =

with ipjtk,P := 2fcd/pV(2fc- - j) the Lp(Rd)-normalized (pre)wavelets, and the

7i,fe,^,p the associated dual functionals. Using the Lp(Rd)-stability of * , we
can replace \\Qk(f)\\pL^Rd) by £ j G Z d |7j,fc,v,p(/)lp ^ obtain

l
'\ (5-8)

with the usual change to a supremum when q = oo. This is the characteri-
zation of the Besov space in terms of wavelet coefficients. When q — p, (5.8)
takes an especially simple form:

fc€Z i/>e* j€Zd

In particular, (5.9) gives an equivalent semi-norm for the Sobolev space

Ha(Rd) by taking p = 2.

6. Applications
6.1. Wavelet compression
We shall present a few examples that indicate how wavelets can be used in
numerical applications. Wavelet techniques have had a particularly signif-
icant impact on data compression. We begin by discussing a problem in
nonlinear approximation that is at the heart of compression algorithms.
 WAVELETS 45

Suppose that * is a (pre)wavelet set and / 6 Lp(Rd), 1 < p < oo, has the
wavelet representation

/ = £ E £ cj
with respect to the Z,p(Rd)-normalized functions ipj,k,P '•= 2fcd/'p^(2fc- — j).
In numerical applications, we must replace the sum in (6.1.1) by a finite
sum, and the question arises as to the most efficient way to accomplish this.
To make this into a well defined mathematical problem, we fix an integer
n, which represents the number of terms we shall allow in the finite sum.
Thus, we want to approximate / in the Lp(M.d) norm by an element from
the set
{ } (6.1.2)
where 4?\fc,V> are arbitrary complex numbers. We have the error of approxi-
mation
an(f)p:= inf \\f - S\\Lp(Rd). (6.1.3)

In contrast to the usual problems in approximation, the set E n is not a linear

space since adding two elements of E n results in an element of £2n> but not
generally an element of S n .
The approximation problem (6.1.3) has a particularly simple solution
when p = 2 and \t is an orthogonal wavelet set. We order the coeffi-
cients Cjifc,v,2(/) by their absolute value. If A n is the set of indices (J, k,xp)
corresponding to n largest values (this set is not necessarily unique), then
Sn = Eo,fc,^)GAn cj>k^,2(f)ipj,k,2 attains the infimum in (6.1.3). For pre-
wavelet sets, this selection is optimal within constants, as follows from the
L2(Kd)-stability of the basis V>j,fc,2-
It is somewhat surprising that the strategy of the preceding paragraph also
is optimal (in a sense to be made clear later) for approximation in Lp(Rd),
p ^ 2. To describe this, wefixa value of 1 < p < oo (slightly weaker results
than those stated later are known when p = oo) and for / 6 Lp(Rd), we let
An denote a set of indices corresponding to n largest values of
We define Sn := E(j,fc»eAn Cj,fc,t/>,p(/)V'j,fc,p and cr n (/) p := ||/ - ^
DeVore et al. (1991d) have established various results that relate <Tn(/)p
with CTn(/)p under certain conditions on ip and the generating function <j>.
For example, it follows from their results that
M / ) p = O(n-<*'d) ^ an(f)p = O(n-a'd) (6.1.4)
for 0 < a < r. Here, the integer r is related to properties of <j>. Namely,
the generating function (j> should satisfy the Strang-Fix conditions (5.4) of
this order and <j> and rp should have sufficient smoothness (for example, Cr
46 R. A. DEVORE AND B. J. LUCIER

is enough). It is also necessary to assume decay for the functions <f> and ip;
sufficiently fast polynomial decay is enough. We caution the reader that the
results in DeVore et al. (1991d) are formulated for one wavelet V and not a
wavelet set \P. However, the same proofs apply in the more general setting.
It is also of interest to characterize the functions / that satisfy (6.1.4).
That is, we would like to know when we can expect the order of approxima-
tion (6.1.4). This has not been accomplished in exactly the form of (6.1.4),
but the following variant has been shown in DeVore et al. (1991d). The
following are equivalent for r := r(a,p) := (a/d+ 1/p)" 1 :

v n (/) P ri < oo,

n=l

f
(ii) f> a / d a n (/) P ] T i < oo,
n=l
(iii) / e B ? ( t T ( R d ) ) . (6.1.5)

Some explanation regarding (6.1.5) is in order. First, ((6.1.5)(i)) is very

close to the condition in (6.1.4). For example, ((6.1.5)(i)) implies (6.1.4),
and if (6.1.4) holds for some /3 > a then ((6.1.5)(i)) is valid. So, roughly
speaking, it is the functions / 6 B"(LT(R(1)) for which the order of approx-
imation in (6.1.4) holds. Second, the characterization (6.1.5) says that it is
those functions with smoothness of order a in XT(Rd), r = (a/d + l/p)~1,
that are approximated with order O(n~a/d) in L p (R d ). This should be con-
trasted with the usual results for approximation from linear spaces (such
as finite element methods), which characterize functions with this approxi-
mation order as having smoothness of order a in Lp(M.d). Since r < p, the
nonlinear approximation problem (6.1.3) provides the approximation order
(6.1.4) for functions with less smoothness than required by linear methods.
The fact that functions with less smoothness can be approximated well by
(6.1.3) is at the essence of wavelet compression. This means that functions
with singularities can be handled numerically. Intuitively this is accom-
plished by retaining in the sum for Sn terms corresponding to functions V'[Link].p
that make a large contribution to / near the singularity. Here the situation
is similar to adaptive methods for piecewise-polynomial (finite element) ap-
proximation that have refined triangulations near a singularity. However,
we want to stress that in wavelet compression, it is simple (almost trivial) to
approximate optimally without encountering problems of triangulation. An
overview of this approach to data compression using wavelets can be found
in DeVore et al. (1991a).
 WAVELETS 47

6.2. Image compression

We explain next how a typical algorithm for compression is implemented
from the theoretical results of the previous section. This has been accom-
plished for surface compression in DeVore et al. (1991b) and image com-
pression in DeVore et al. (1991c). We shall discuss only image compression.
A digitized grey-scale image consists of an array of picture elements (pix-
els) represented by numbers that correspond to the brightness (grey scale)
of each pixel, with 0 = black, say, and 255 = white. A grey-scale image
has, say, 1024 x 1024 such numbers taking integer values between 0 and 255.
Thus, the image is given by a matrix 0>j)je{o,...,iO23}2 with Pj e {0, • • •, 255}.
It would take a data file of about one million bytes to encode such an image.
For purposes of transmission or storage, it is desirable to compress this file.
To use wavelets for image compression, we proceed as follows. We think
of the pixel values as associated with the points j2~m, j e [0,2 m ) 2 , m = 10,
of the unit square [0,1]2. In this way, we can think of the image as a
discretization of a function / defined on this square.
We choose a function <t> satisfying the assumptions of multi-resolution, and
a corresponding wavelet set W that provides a stable basis for L2OR2). Thus,
^ would consist of three functions, which we shall assume are of compact
support.
We would like to represent the image as a wavelet sum. For this purpose,
we select coefficients 7j and consider the function

J60

with ft the set of indices for which <pj>m does not vanish identically on [0,1]2.
We think of / as the image and apply the results of the preceding section
to compress / .
The coefficients fj are to be determined numerically from the pixel values;
choosing good values of jj is a nontrivial problem, which we do not discuss.
A typical choice is to take fj = pj for j2~m G [0,1]2 and some extension of
these values for other j .
Once the coefficients (jj) have been determined, we use a fast wavelet
transform to write / in its wavelet decomposition

(6-2-2)

with respect to the Z2(M2)-normalized ipj^ '•= 2kip(2k-—j). We can find the
coefficients of / with respect to the Lp(R2)-normalized V's by the relation
c
j,k,ip,p = 22fc(1/p~1/2)cJ)fc)y,. The projection Pof has very few terms, which
we take intact into the compressed representation at little cost.
To apply the compression algorithm of the previous section, we need to
48 R. A. DEVORE AND B. J. LUCIER

decide on a suitable norm in which to measure the error. The Z<2(K2) norm is
most commonly used, but we argue in DeVore et al. (1991c) that the Li(R 2 )
is a better model for the human eye for error with high spatial frequency.
If one decides to use the LP(R2) norm to measure compression error, then
the algorithm of the previous section orders the Lp(ffi2)-normalized wavelet
coefficients Cjtk,ip,P and chooses the largest of these to keep. Optimally, one
would send coefficients in decreasing order of size. Thus, we find a (small)
set A of ordered triples {(j,k,ip)} that index the largest values of
and use for our compressed image

This sum has |A| terms.

This method of sending coefficients sequentially across a communications
link to allow gradual reconstruction of an image by the receiver is known
as progressive transmission. Our criterion provides a new ordering for the
coefficients to be transmitted that depends on the value of p. However,
sorting the coefficients requires O(m22m) operations, while the fast wavelet
transform itself takes but O(22m) operations. Thus, a faster compression
method that does not rely on sorting is to be preferred; we proceed to give
one.
We discuss compression in L2OR2) for a moment. As noted earlier, the
optimal algorithm is to keep the largest Z,2(R2) coefficients and to discard
the other coefficients. The coefficients to keep can be determined by fixing
any the following quantities:
(i) JV:=|A|,
/ \l/2
2
(ii) 11/ - gWw = { £ )
v
UM)
(iii) e:= inf \cjM\.
Ofc^)€A
Setting any one of these quantities determines the other two, and by exten-
sion the set A, for any function / . In other words, we can prescribe either
the number of nonzero coefficients N, the total error ||/ — <7||L2(R2)> or e,
which we consider to be a measure of the local error. If we determine which
triples (j, k, ip) to include in A by the third criterion, then we do not need
to sort the coefficients, for we can sequentially examine each coefficient and
put (j,k,ip) into A whenever I c ^ ^ l > e. This is known as threshold coding
to the engineers, because one keeps only those coefficients that exceed a
specified threshold.
Even more compression can be achieved by noting that we should keep
only the most significant bits of the coefficients Cjtk,<j>,p- Thus, we choose a
 WAVELETS 49

tolerance e > 0 and we take in the compressed approximation for each ipj,k,p
a coefficient Cjtkj>,p such that
\cj,k,yl>,p - cj,k,i>,p\ < e
(6.2.3)
with the proviso that Cj,fc,^,p — 0 whenever \CJ^,P\ < e. Then S :=
Sj,fc,^ Cj^xijffiipj^p represents our compressed image, and (6.1.5) holds for
this approximation.
This process of keeping only the most significant bits of Cjtk,4>,P is known
in the engineering literature as scalar quantization. The dependence on the
dyadic level k and the space Lp([0,I]2) in which the error is measured is
brought out more clearly when using £00(M2)-normalized wavelets, i.e. when
i>j,k — V'(2fc- — j)- For these wavelets, as the dyadic level increases, the
number of bits of c,^^ taken in Cj^,^ decreases. For example, if p = 2, we
would take one less bit at each increment of the dyadic level. On the other
hand, if the compression is done in the L\([0,1]2) norm, than we would take
two fewer bits as we change dyadic levels. See DeVore et al. (1991c).

6.3. The numerical solution of partial differential equations

Wavelets are currently being investigated for the numerical solution of dif-
ferential and integral equations (see, e.g., the papers of Beylkin et al. (1991)
and Jaffard (1991)). While these applications are only now being developed,
we shall consider a couple of simple examples to illustrate the potential of
wavelets in this direction.
Elliptic equations The Galerkin method applied to elliptic partial dif-
ferential equations gives rise to a matrix problem that involves a so-called
stiffness matrix. We present a simple example that illustrates the perhaps
surprising fact that the stiffness matrix derived from the wavelet basis can
be preconditioned trivially to have a uniformly bounded condition number.
In general, this property allows one to use iterative methods, such as the
conjugate gradient method, to solve linear systems with great efficiency. The
linear systems that arise by discretizing elliptic PDEs have a lot of structure
and can in no way be considered general linear systems, and there are many
very efficient numerical methods, such as multi-grid, that exploit the special
structure of these linear systems to solve these systems to high accuracy
with very low operation counts. We do not yet know of a complete analysis
that shows that computations with wavelets can be more efficient than ex-
isting multi-grid methods when applied to the linear systems that arise by
discretizing elliptic PDEs in the usual way.
Rather than consider Dirichlet and Neumann boundary value problems
in several space dimensions, as discussed in Jaffard (1991), we shall present
only a simple, periodic, second-order ODE that illustrates the main points.
We shall consider functions defined on the one-dimensional torus T, which
50 R. A. DEVORE AND B. J. LUCIER

is equivalent to [0,1] with the end-points identified, and search for u = u(x),
x € T, that satisfies the equation
- u"(x) + u(x) = f(x), x G T, (6.3.1)
l
with / € L2(T). In variational form, the solution u € W (L2(T)) of (6.3.1)
satisfies
/ (u'v' + uv)= f fv, (6.3.2)
h JT
for all v € W1(L2(T)). We remark that we can take

To approximate u by Galerkin's method, we must choose a finite-dimen-

sional subspace of W1(L,2(T)), which we shall choose to be a space spanned
by wavelets defined on the circle T. We indicate briefly how to construct
periodic wavelets on T from wavelets on R.
For any / with compact support defined on R, the function

is a function of period 1, which we call the periodization of / . To obtain

wavelets on T, we apply this periodization to wavelets on R. To be specific,
we consider only the Daubechies wavelets ip := T>2N with N > 2 because
they are contained in W1(L,2(T)). (The first nontrivial Daubechies wavelet,
V4, is in Wa(L2(T)) for all a < 1 (see Eirola (1992)), but we do not know if
it is in W1(L,2(T)).) Let 4> be the function of Section 3.5 that gives rise to ip.
For each k € Z, we let <j>°jk and ipjk be the periodization of the functions 0 ^
and ipj^ respectively. We define Sk to be the linear span of the functions
tf>° , j = 0 , . . . , 2* — 1. The functions in this space are clearly of period
one. We also note that the <j>jk are orthogonal. Indeed, because they are
periodic,

= / <t>jjt4>°i',k = E / ^j,k<t>j',k( • +t) = E / (t>j,k(t>j'+e2l<,k-

eez
If j ^ j ' , each integral in the last sum is zero, since we never have j =
j ' + 2k£. If j = f, then exactly one integral in the last sum is nonzero and
its value is one. Similarly, we find that ipjk, j = 0 , . . . , 2k — 1, k > 0, is an
orthonormal system for Z^T). It is easy to check that by adjoining 0o O ,
which is identically one on T, this orthonormal system is complete.
Returning to our construction of periodic wavelet spaces, we define Wk
 WAVELETS 51

to be the linear span of the functions V'j,*) j = 0 , . . . , 2* — 1. Then Sk+1 =

Sk © Wk. To simplify our notation in this section, in which we refer to
periodic wavelets only, we drop the superscripts and tildes, and denote by
tf>j;/(., i/tjjc, Sk, and Wk, the periodic wavelet bases and spaces.
Returning to the numerical solution of (6.3.1), we choose a positive value
of m and approximate u by an element um € Sm that satisfies:

/ « y + umv) = f fv, v£Sm. (6.3.3)

JT JT
We can write um and v as in (6.2.2). For example,
m-l m-12fc-l
«m = Poum + E ^ U m =
"T^o.° + E L W ' fc
)^'.fc- (6.3.4)
*=0 fe=0 j=0
Because <^o,o is identically one on the circle T, 7 is just the average of um
on T. From (6.3.3) we see that with v = 1, 7 = / T / . Thus, we need to
determine /3(j, fc).
If we replace um in (6.3.3) by its representation (6.3.4), we arrive at a
system of equations
TO_12*-i

EE
fc=0 j=0

for f = 0 , . . . , 2k> — 1 and k' = 0 , . . . , m — 1, or, more succinctly,

T(3 = f (6.3.5)
where the typical entry in T is

li
11
and / a vector with components

J jW, P = G0(j'i *))j=0,...,2*-l, fc=0,...,m-l

is the coemcient vector of the unknown function.

The convergence rate of the conjugate gradient method, a popular itera-
tive method for the solution of systems like (6.3.5), depends on the condi-
tion number, K(T) = j|y || HT1"11|, of the symmetric, positive definite, stiffness
matrix T. Now,
||T||= sup <x*Ta and WT^W'1 = mi ct*Ta,
||a|| < 2 =l l|a|l«3=i

where a := (a(j,fc))j=o,...,2fc-i,fc=o,...,m-i-For any vector a , we form the

52 R. A. DEVORE AND B. J. LUCIER

function S := S(a) 6 <Sm by

m-12*-l

fe=O j=0

For example, um = 7<fo,0 + S(P)- It follows easily that

fl m-12k-l
2
\a*T<*\ = I ([Sf + S ) = ||5||^ 1 ( L 2 ( R d ) ) « £ £ |2*a(j,*)l2- (6-3.6)

The last equivalence in (6.3.6) is a variant of (5.9) (for q = p = 2). In

(5.9), we considered wavelet representations beginning at the dyadic level
k = — oo; we could have just as easily begin at the level k = 0 (by including
) and arrived at the equivalence in (6.3.6).
Equation (6.3.6) shows that the matrix DTD, where D has entries

satisfies
m-12*-l

i.e. the stiffness matrix T can be trivially preconditioned to have a condition

number K(DTD) := ||J9TjD||||.D~ 1 r~ 1 .D- 1 || that is bounded independently
of m.
Because the subspace Sm of 1/2(1) generated by the Daubechies wavelet
£>2jv locally contains all polynomials of degree < N, we have that

Because of the special form of (6.3.1), um is in fact the W1(L2(T)) projection

of u onto Sm, so we have

One can treat in an almost identical way the equation (6.3.1) defined on
the torus Trf with the left-hand side replaced by — V • (aVu) + bu with a and
6 bounded, smooth, positive functions on T d , and obtain the same uniform
bound on the condition number. Of course, ultimately, one would like to
handle elliptic boundary value problems for a domain SI. First results in
this direction have been obtained by Jaffard (1991) and in a slightly differ-
ent (nonwavelet) setting by Oswald (see for example Oswald (1991)). For
example, Jaffard's approach to elliptic equations with Dirichlet boundary
conditions is to transform the equation to one with zero boundary condi-
tions by extending the boundary function into the interior of the domain Q.
He then employs in a Galerkin method wavelets whose support is contained
 WAVELETS 53

strictly inside fi. However, there has not yet been an analysis of the desired
relationship between the extension and the wavelets employed. Another ap-
proach is to develop wavelets for the given domain (which do not vanish on
the boundary) (see Jaifard and Meyer (1989)).
Employing wavelets for elliptic problems as outlined earlier is similar to
the use of hierarchical bases in the context of multi-grid. However, two
points suggest that wavelet bases may be more useful than hierarchical bases.
First, the wavelet bases are L2(Kd)-stable, while the hierarchical basis are
not. Second, one can choose from a much greater variety of wavelet bases
with various approximation properties (values of N).
Finally, we mention the great potential for compression to be used in
conjunction with wavelets for elliptic problems in a similar way to adaptive
finite elements (see also Jaffard (1991)).

6.4- Time-dependent problems

Wavelets have a potential application for the numerical solution of time-
dependent parabolic and hyperbolic problems. We mention one particular
application where the potential of wavelets has at least a theoretical foun-
dation.
We consider the solution u{x, t) of the scalar hyperbolic conservation law
ut + f{u)x = 0, x e R, t> 0,
u(x,0) = uo(x), x€R, (6.4.1)
in one space dimension. It is well known that the solution to (6.4.1) develops
discontinuities (called 'shocks') even when the initial condition is smooth.
This makes the numerical treatment of (6.4.1), and even more so its analogue
in several space dimensions, somewhat subtle. The appearance of shocks
calls for adaptive or nonlinear methods.
Considering the appearance of discontinuities in the solution u to (6.4.1),
the following regularity result of DeVore and Lucier (1990, 1988), is quite
surprising. If the flux / in (6.4.1) is strictly convex and suitably smooth,
and UQ has bounded variation, then it has been shown that, for any a > 0
and T := r(a) := (a + I ) " 1 ,
u0 € B?(Lr(R)) => u( •, t) G B?(LT(R)), (6.4.2)
for all later time t > 0. That is, if UQ has smoothness of order a in LT(R)
then so will u( •, t) for all later time t > 0.
The regularity result (6.4.2) has an interpretation in terms of wavelet de-
compositions. We cannot use orthogonal wavelets in these decompositions
because, as we pointed out earlier, they do not provide stable representations
of functions in L\ (R) because they have mean value zero. However, the char-
acterization of Besov spaces can be carried out using other, nonorthogonal,
54 R. A. D E V O R E AND B. J. LUCIER

wavelets, such as B-splines. With this caveat, (6.4.2) says that whenever
uo has a wavelet decomposition X^[Link]^o)^,*;,!: with certain Li(R)-
normalized wavelets ip, whose coefficients satisfy
<°o, (6.4.3)
fcez jez
then u( •, t) has a similar wavelet decomposition with the same control
(6.4.3) on the coefficients. We want to stress that the results in DeVore
and Lucier (1990,1988) do not describe directly how to determine wavelet
coefficients at a later time t > 0 from those of the initial function «o- That
is, there is no direct, theoretically correct, numerical method known to us
that describes how to update coefficients with time so that (6.4.3) holds.
The regularity result (6.4.2) is proved by showing that whenever UQ can be
approximated well in Li(R) by piecewise polynomials with free (variable)
knots, then u(- ,t) can be approximated in the same norm by piecewise
algebraic functions (of a certain type) with free knots.
Finally, we mention that the authors have also shown that the analogue
of regularity result (6.4.3) does not hold in more than one space dimension.
From the point of view of wavelet decompositions, our results seem to indi-
cate that the wavelets described in this presentation are too symmetric to
effectively handle the diverse types of singularities that arise in the solution
of conservation laws in several space dimensions.

REFERENCES
G. Battle (1987), 'A block spin construction of ondelettes, Part I: Lemarie func-
tions', Comm. Math. Phys. 110, 601-615.
A. Ben Artzi and A. Ron (1990), 'On the integer translates of a compactly supported
function: Dual bases and linear projectors', SIAM J. Math. Anal. 21, 1550-
1562.
G. Beylkin, R. Coifman and V. Rokhlin (1991), 'Fast wavelet transforms and nu-
merical algorithms I', Comm. Pure Appl. Math. XLIV, 141-183.
C. de Boor, R. DeVore and A. Ron (1991a), 'Approximation from shift invariant
spaces', Preprint.
C. de Boor, R. DeVore and A. Ron (1991b), 'Wavelets and pre-wavelets', Preprint.
C. de Boor, R. DeVore and A. Ron (1991c), 'The structure of finitely generated
shift invariant spaces in Z,2(Md)', Preprint.
C. de Boor and A. Ron (1991), 'Fourier analysis of approximation orders from
principal shift invariant spaces', Constructive Approx., to appear.
A. Cavaretta, W. Dahmen and C. Micchelli (1991), 'Subdivision for computer aided
geometric design', Memoirs Amer. Math. Soc. 93.
C. K. Chui, J. Stockier and J. D. Ward (1991), 'Compactly supported box-spline
wavelets', CAT Report, Texas A&M University 230
C. K. Chui and J. Z. Wang (1990), 'A general framework for compactly supported
splines and wavelets', CAT Report, Texas A&M University 219.
 WAVELETS 55

C. K. Chui and J. Z. Wang (1991), 'On compactly supported spline wavelets and a
duality principle', Trans. Amer. Math. Soc, to appear.
W. Dahmen and C. Micchelli (1990), 'On stationary subdivision and the construc-
tion of compactly supported wavelets', Multivariate Approximation and Inter-
polation Vol. 94 (W. Haussmann and K. Jetter, eds) ISNM, Birkhauser Verlag
(Boston) pp. 69-89.
I. Daubechies (1988), 'Orthonormal bases of compactly supported wavelets', Comm.
Pure Appl. Math. XLI, 909-996.
I. Daubechies (1990), 'The wavelet transform, time-frequency localization and sig-
nal analysis', IEEE Trans. Information Theory 36, 961-1005.
R. DeVore, B. Jawerth and B. Lucier (1991a), 'Data compression using wavelets:
Error, smoothness, and quantization', DCC-91, Data Compression Confer-
ence (A. Storer and J. H. Reif, eds), IEEE Computer Society Press (Los
Alamitos, CA) pp. 186-195.
R. DeVore, B. Jawerth and B. Lucier (1991b), 'Surface compression', Computer
Aided Geometric Design, to appear.
R. DeVore, B. Jawerth, and B. Lucier (1991c), 'Image compression through wavelet
transform coding', IEEE Trans. Information Theory to appear.
R. DeVore, B. Jawerth and V. Popov (1991d), 'Compression of wavelet decompo-
sitions', Amer. J. Math., to appear.
R. DeVore and G. Lorentz (1992), Constructive Approximation Springer
Grundlehren (Heidelberg).
R. DeVore and B. Lucier (1989), 'High order regularity for solutions of the inviscid
Burgers equation', Nonlinear Hyperbolic Problems (Proc. Advanced Research
Workshop, Bordeaux, France, June 1988, Springer Lecture Notes in Math-
ematics, 1402) (C. Carasso, P. Charrier, B. Hanouzet and J.-L. Joly, eds)
Springer-Verlag (New York) pp. 147-154.
R. DeVore and B. Lucier (1990), 'High order regularity for conservation laws',
Indiana Univ. Math. J. 39, 413-430.
R. DeVore and V. Popov (1988), 'Interpolation spaces and non-linear approxima-
tion', Function Spaces and Applications (Springer Lecture Notes in Mathe-
matics, 1302) (M. Cwikel, J. Peetre, Y. Sagher and H. Wallin, eds) Springer
(New York) pp. 191-205.
T. Eirola (1992) 'Sobolev characterization of solutions of dilation equations', SIAM
J. Math. Anal, to appear.
M. Frazier and B. Jawerth (1990), 'A discrete transform and decompositions of
distribution spaces', J. Funct. Anal. 93, 34-170.
M. Frazier, B. Jawerth and G. Weiss (1991), Littlewood-Paley Theory and the Study
of Function Spaces Conference Board of the Mathematical Sciences, Regional
Conference Series in Mathematics, Number 79, Amer. Math. Soc. (Providence,
RI)
S. Jaffard (1991) 'Wavelet methods for fast resolution of elliptic problems', Preprint
S. Jaflfard and Y. Meyer (1989), 'Bases ondelettes dans des ouverts de Mn>, J. Math.
Pures Appl. 68, 95-108.
R-Q. Jia and C. Micchelli (1991), 'Using the refinement equation for the construc-
tion of pre-wavelets II: power of two', Curves and Surfaces, (P. J. Laurent, A.
LeM4haute, and L. Schumaker, eds) Academic Press (New York) pp. 209-246.
 56 R. A. D E V O R E AND B. J. LUCIER

S. Karlin and W. Studden (1966), Tchebycheff Systems with Application to Analysis

and Statistics Wiley Interscience (New York).
R. A. H. Lorentz and W. R. Madych (1991) Wavelets and generalized box splines,
preprint.
S. G. Mallat (1989), 'Multi-resolution approximations and wavelet orthonormal
bases of £ 2 ( R ) \ Trans. Amer. Math. Soc. 315, 69-87.
Y. Meyer (1987), 'Wavelets with compact support', Zygmund Lectures (University
of Chicago).
Y. Meyer (1990), Ondelettes et Operateurs I: Ondelettes Hermann (Paris).
C. Micchelli (1991), 'Using the refinement equation for the construction of pre-
wavelets', Numer. Alg. 1, 75-116.
P. Oswald (1991) 'On discrete norm estimates related to multi-level preconditioners
in the finite elment method', Preprint.
S. D. Riemenschneider and Z. W. Shen (1991), 'Box splines, cardinal splines, and
wavelets', Approximation Theory and Functional Analysis (C. K. Chui, ed.)
Academic Press (New York) pp. 133-149.
S. D. Riemenschneider and Z. W. Shen (1992) 'Wavelets and pre-wavelets in low
dimensions', J. Approx. Theory to appear.
0 . Rioul and M. Vetterli (1991), 'Wavelets and signal processing', IEEE Signal
Processing Magazine 8, (4) 14-38.
1. Schoenberg (1946), 'Contributions to the problem of approximation of equidistant
data by analytic functions, Parts A & B', Quart. Appl. Math. IV, 45-99,112-
141.
I. Schoenberg (1973), Cardinal Spline Interpolation (Regional Conference Series in
Applied Mathematics 12) SIAM (Philadelphia).
G. Strang and G. Fix (1973), 'A Fourier analysis of the finite element variations!
method', C.I.M.E. II Ciclo 1971, Constructive Aspects of Functional Analysis
(G. Geymonat, ed.) pp. 793-840.
J. O. Stromberg (1981) 'A modified Franklin system and higher-order spline systems
on K n as unconditional basis for Hardy spaces', Conference in Harmonic Anal-
ysis in Honor of Antoni Zygmund Vol II (W. Beckner et al., eds) Wadsworth
pp. 475-493.

View publication stats