Oblique and Hierarchical Multiwavelet Bases

APPLIED AND COMPUTATIONAL HARMONIC ANALYSIS ARTICLE NO.

4, 231–263 (1997)

HA970211

Oblique and Hierarchical Multiwavelet Bases AKRAM ALDROUBI1 Vanderbilt University, Department of Mathematics, 1326 Stevenson Center, Nashville, Tennessee 37240-0001

Communicated by Pierre G. Lemarie´-Rieusset Received March 26, 1996; revised December 30, 1996

We develop the theory of oblique multiwavelet bases, which encompasses the orthogonal, semiorthogonal, and biorthogonal cases, and we circumvent the noncommutativity problems that arise in the construction of multiwavelets. Oblique multiwavelets preserve the advantages of orthogonal and biorthogonal wavelets and enhance the flexibility of the theory to accommodate a wider variety of wavelet bases. For example, for a given multiresolution, we can construct supercompact wavelets for which the support is half the size of the shortest orthogonal, semiorthogonal, or biorthogonal wavelet. The theory also produces the h-type, piecewise linear hierarchical bases used in finite element methods, and it allows us to construct new h-type, smooth hierarchical bases, as well as h-type hierarchical bases that use several template functions. For the hierarchical bases, and for all other types of oblique wavelets, the expansion of a function can still be implemented with a perfect reconstruction filter bank. We illustrate the results using the Haar scaling function and the Cohen–Daubechies–Plonka multiscaling function. We also construct a supercompact spline uniwavelet of order 3 and a hierarchical basis that is based on the Hermit cubic spline, and we explicitly give the coefficients of the corresponding filter bank. q 1997 Academic Press

1. INTRODUCTION

A multiscaling function is a vector function F(x) Å (f1(x), f2(x), f3(x), . . . , fr(x))T that can be used to generate the multiresolution spaces

H H∑ r

Vj Å

∑ ∑ cij(k)20j/2fi

iÅ1 k√Z

Å

S D

x 0 k ; cij √ l2 , i Å 1, . . . , r 2j

CjT(k)Fj(x 0 2jk); Cj(k) √ lr2 ,

k√Z

1

J

J (1.1)

E-mail: [email protected].

231 1063-5203/97 $25.00 Copyright q 1997 by Academic Press All rights of reproduction in any form reserved.

6111$$0211

06-10-97 11:03:28

achaa

AP: ACHA

232

AKRAM ALDROUBI

where Cj(k) is the vector Cj(k) Å (c1j (k), c2j (k), . . . , crj(k))T, T denotes the matrix transpose operator, lr2 Å l2 1 rrr 1 l2 , and the vector Fj(x) is defined to be Fj(x) Å (f1j (x), f2j (x), . . . , frj(x))T

S SD SD

Å 20j/2 f1

S DD

x x x , f2 j , . . . , fr j 2j 2 2

T

.

(1.2)

The spaces Vj must satisfy all the properties of multiresolutions [24]. In particular, they must be closed, nested: rrrV2 , V1 , V0 , V01 , V02 , rrr,

and they must satisfy the two size properties V0` Å < Vj Å L2 , j√Z

V` Å > Vj Å {0}. j√Z

Moreover, the set {fij,k(x) Å fij(x 0 2jk); k √ Z, i Å 1, . . . , r.}, must form a Riesz basis of Vj . From the nestedness property, it follows that the multiscaling function F must satisfy the vector two-scale equation F

SD

x Å 2 ∑ H1(k)F(x 0 k), 2 k√Z

(1.3)

where H1(k) is an r 1 r matrix-sequence called the generating sequence. It is often desirable to construct scaling functions by starting from a generating sequence H1(k). For this purpose, we can take the Fourier transform of (1.3) and obtain the Fourier relation FO (2f ) Å HO 1( f )FO ( f ),

(1.4)

where FO ( f ) is the Fourier transform of the components of F: [FO ( f )]j Å fO j( f ) Å

* f (t)e j

0i2pft

dt,

R

ˆ 1( f ) is the Fourier transform of H1(k). It is given by and where H HO 1( f ) Å ∑ H1(k)e0i2pfk. k√Z

6111$$0211

06-10-97 11:03:28

achaa

AP: ACHA

233

OBLIQUE MULTIWAVELETS

The repeated application of the relation in (1.4) yields

SD S D

FO ( f ) Å HO 1

SDSD

f f f f HO 1 2 rrrHO 1 j FO j . 2 2 2 2

This last equality suggests a construction of F in Fourier domain. The candidate function FO ( f ) is an infinite product of the form

SD S D

FO ( f ) Å lim HO 1 jr`

SD

f f f HO 1 2 rrrHO 1 j a, 2 2 2

(1.5)

where a is a vector in Rr. For the infinite product in (1.5) to converge to a continuous ˆ 1(0) must have the eigenvalue l Å 1, and a must be an associated right function, H eigenvector [9, 18, 19]. There are other conditions that must be met, as well. A ˆ 1(0) ˆ 1( f ) 0 H sufficient condition is to have an upper-bound on the matrix-norm of H of the form \HO 1( f ) 0 HO 1(0)\ £ CÉfÉa,

ˆ 1(0) of for some constant a ú 0, and an upper-bound on the spectral radius r0 of H a the form r0 £ 2 [9]. Other sufficient conditions can be found in [18, 20]. The convergence of the infinite product in (1.5) is still not sufficient to ensure that it produces a scaling function. Other conditions must be satisfied. In particular, the conditions on H1(k) that coerce V0` to be dense in L2 with approximation order m have been discussed by Heil, Strela, and Strang as well as by Cohen, Daubechies, and Plonka [9, 19, 28]. Finally, for {fij,k; k √ Z, i Å 1, . . . , r} to form a Riesz basis of Vj , an additional condition must be satisfied on the autocorrelation matrix-sequence [A]m,n(k) Å »fm(x), fn(x 0 k)…L2 .

(1.6)

Specifically, the set {fij,k; k √ Z, i Å 1, . . . , r} forms a Riesz basis, if and only if, the Fourier transform AO ( f ) Å ∑ FO ( f / k)FO *( f / k)

(1.7)

k√Z

of A(k) is positive definite (‘‘∗’’ in (1.7) denotes the Hermitian transpose operator), self-adjoint, and there exist two positive constants m ú 0 and M ú 0 such that ∀Y √ Rr: m\Y\2Rr £ »AO ( f )Y, Y…L2 £ M\Y\2Rr

a.e. f.

(1.8)

This result can be found in [12, 13, 16], and it is a special case of the result in [1, Theorem 2.1].

6111$$0211

06-10-97 11:03:28

achaa

AP: ACHA

234

AKRAM ALDROUBI

1.1. Wavelet Spaces and Multiwavelets By definition, a sequence of spaces {Wj}j√Z are called wavelet spaces associated with the multiresolution spaces if and only if • The spaces {Wj}j√Z complement the spaces {Vj}j√Z :

Vj/1 / Wj/1 Å Vj .

(1.9)

• The spaces {Wj}j√Z are generated by shifts and dilations of a vector function C(x) Å (c1(x), c2(x), . . . , cr(x))T called a multiwavelet. • The set

{cij,k(x) Å 20j/2c1j (x 0 2jk), . . . , 20j/2crj(x 0 2jk); k √ Z} forms a Riesz basis of

H∑ ∑ r

Wj Å

dij(k)20j/2ci

iÅ1 k√Z

S D

x 0 k ; dij √ l2, i Å 1, . . . , r. 2j

J

Å { ∑ DjT(k)Cj(x 0 2jk); Dj(k) √ lr2},

(1.10)

k√Z

where Cj(x) is defined to be Cj(x) :Å 20j/2C(x/2j).

(1.11)

Note that we do not require orthogonality between the spaces Vj and Wj . Therefore, the spaces Wj are not necessarily orthogonal to each other, either. Requirement (1.9) implies that the multiwavelet must be a linear combination of shifts and dilations of the scaling function. Therefore, we are led to the relation C

SD

x Å 2 ∑ (d1∗G1)(k)F(x 0 k), 2 k√Z

(1.12)

where G1(k) is a matrix-sequence, di(k) is the unit impulse sequence located at k Å i,

di(k) Å

H

1, k Å i

(1.13)

0, k x i;

and the generalized convolution (a∗B)(k) between a scalar-sequence a(k) and a matrixsequence B(k) is the matrix-sequence C(k) defined by Ci,j(k) Å ∑ a(l)Bi,j(k 0 l). l√Z

6111$$0211

06-10-97 11:03:28

achaa

AP: ACHA

235

OBLIQUE MULTIWAVELETS

From this definition, it follows that the unit impulse d1(k) in (1.12) is used to shift the sequence G1(k), i.e., (d1∗G1)(k) Å G1(k 0 1). In the context of this paper, we will also use the symbol ‘‘∗’’ to denote another generalized convolution between two sequences of matrices. Specifically, the convolution C(k) Å (A∗B)(k) between the m 1 r matrixsequence {A(k)}k√Z and the r 1 n matrix-sequence {B(k)}k√Z is the m 1 n matrixsequence, defined in terms of the convolution between the entries of A and B as r

Ci,j(k) :Å ∑ ∑ Ai,l(h)Bl,j(k 0 h).

(1.14)

lÅ1 h√Z

Remark 1.1. The generalized convolution (1.14) is a combination of regular convolution of sequences and matrix multiplication. In particular, if m Å r Å n Å 1, then all the matrices become scalar sequences, and the generalized convolution reduces to the usual convolution between sequences. 1.2. Classification of Wavelet Bases According to the angles between the spaces {Wj}j√Z and the choice of bases for these spaces, wavelets, and multiwavelets can be classified into several categories: • Orthogonal wavelet and multiwavelet bases [11, 13, 23, 24, 31]. For this case, there are two conditions that must be satisfied:

(1) The wavelet spaces must be orthogonal to each other: W j ⊥ Wl

∀j x l.

(2) The set {cij,k ; k √ Z, i Å 1, . . . , r.} must be an orthogonal basis of Wj . • Semiorthogonal wavelet and multiwavelet bases [3, 6, 12, 27, 29, 33]. For these bases the only requirement is that Wj ⊥ Wl ∀j x l. In this case, {cij,k ; k √ Z, i Å 1, . . . , r.} is not necessarily an orthogonal basis of Wj . • biorthogonal wavelet bases [8, 25]. This case consists of a pair of wavelet bases {cij,k ; k √ Z, i Å 1, . . . , r.} and {cH ij,k ; k √ Z, i Å 1, . . . , r.} generating a pair of ˜ j and satisfying the biorthogonality condition wavelet spaces Wj and W »cij,k , cH lm,n…L2 Å d0(j 0 m)d0(k 0 n)d0(i 0 l).

(1.15)

• Oblique wavelet and multiwavelet bases [2]. In this case, we only require the set {cij,k ; k √ Z, i Å 1, . . . , r.} to be a Riesz basis of Wj . We do not require orthogonality between the wavelet spaces, as in the orthogonal and semiorthogonal cases; and unlike the biorthogonal case, we do not require associated wavelet spaces ˜ j or the existence of a biorthogonal wavelet. W

For the case of uniwavelets, a good overview can be found in [21]. The first orthogonal multiwavelet bases were introduced by Donovan, Geronimo, Hardin, and Massopust [13]. These bases have the remarkable properties of being orthogonal, regular, compactly supported, and symmetrical. It was not possible to construct wavelet

6111$$0211

06-10-97 11:03:28

achaa

AP: ACHA

236

AKRAM ALDROUBI

bases (r Å 1) satisfying these four properties before [11]. Other orthogonal multiwavelet bases have also been found [31]. Semiorthogonal constructions of multiwavelets were first introduced in [15, 16], and more recent work can be found in [29]. Oblique multiwavelet bases were introduced in [2], and the underlying theory encompasses the theory of biorthogonal multiwavelets of which recent examples have been found in [25]. The theory of oblique multiwavelets also produces the piecewise linear, htype hierarchical basis used in finite element methods [4, 38], and it allows us to extend the construction to new, smoother hierarchical bases as well as to the hierarchical bases that are based on multiple template functions. What is remarkable is that the expansion algorithm in terms of these bases can also be obtained by filter banks. As an example, in Section 7.3, we construct a hierarchical basis that is based on the Hermite cubic splines, and we give the coefficients of the associated filter bank explicitly. In this paper, we build a multiwavelet theory in which the wavelet spaces are merely complementary to the multiresolution spaces, resulting in the oblique multiwavelets. Moreover, the computation of the oblique multiwavelet transform can be calculated using fast filter-bank algorithms with a complexity of O(N) for a sequence of length N. This construction allows us more freedom in choosing the properties of the wavelet spaces. At the same time, we have enough control on the relationship between the multiresolution spaces and the wavelet spaces to provide fast and stable filter-bank algorithms. The paper is organized as follows: In Section 2, we present and discuss a simple example based on the Haar multiresolution. Section 3 is devoted to the development of the theory of oblique multiwavelet bases. The special case of biorthogonal multiwavelets is discussed in Section 3.2. The implementation of the multiwavelet transform is discussed in Section 4, where we also introduce the concept of perfect demultiplexing filter banks. In Section 5, we discuss the notion of equivalence between wavelet bases. In Section 6, we construct several examples of oblique wavelet and multiwavelet bases. In particular, we use the multiscaling function of Cohen, Daubechies, and Plonka to construct several examples. Finally, in Section 7, we discuss the connection with hierarchical bases and construct a supercompact hierarchical wavelet basis and a Hermite cubic splines hierarchical basis that uses two template functions. 2. A HAAR OBLIQUE WAVELET

Consider a function g which belongs to the space of piecewise constant functions V0 Å {( c0(k)x[0,1)(x 0 k); c0 √ l2}, g(x) Å ∑ c0(k)x[0,1)(x 0 k),

(2.1)

k√Z

where x[a,b) is the characteristic function on the interval [a, b) (i.e., x[a,b)(x) Å 1 ∀x √ [a, b), and x[a,b)(x) Å 0 ∀x √ / [a, b). Using the fact that x[0,2) Å x[0,1) / x[1,2) , we see that g can be written as g(x) Å ∑ c0(2k)x[0,2)(x 0 2k) / ∑ (c0(2k / 1) 0 c0(2k))x[1,2)(x 0 2k). k√Z

6111$$0211

k√Z

06-10-97 11:03:28

achaa

AP: ACHA

(2.2)

237

OBLIQUE MULTIWAVELETS

Thus, we have decomposed an arbitrary piecewise constant function g √ V0 in terms of two functions, gV1(x) Å ∑ c0(2k)x[0,2)(x 0 2k) and gW1(x) Å ∑ (c0(2k / 1) 0 c0(2k))x[1,2)(x 0 2k). The function gV1 belongs to the piecewise constant space V1 generated by the Riesz basis {f1,k(x) Å 201/2x[0,2)(x 0 2k); k √ Z}, while the function gW1 belongs to the space W1 , defined by W1 Å { ∑ d1(k)x[1,2)(x 0 2k); d1 √ l2}.

(2.3)

k√Z

It is not difficult to prove that the set {c1,k(x) Å x[1,2)(x 0 2k); k √ Z} is a Riesz basis of W1 . Before we proceed further, we wish to highlight the fact that any function w1 in W1 must vanish on the intervals [2k, 2k / 1) for all k √ Z. This property is due to the fact that the shifts of x[1,2) in (2.3) are multiples of 2 while the support of x[1,2) has length 1. Clearly, V1 , V0 and W1 , V0 . Thus, we can apply the same decomposition algorithm to gV1 and gW1 , now both viewed as elements of V0 . By inspection, we see that this procedure leaves gV1 and gW1 invariant. Moreover, we have that the intersection V1 > W1 Å {0}. To see this, we simply note that if w1 √ W1 , then w1(2k / 12) Å 0 for any k √ Z. On the other hand, if w1(x) Å (l c1(l)x[0,2)(x 0 2l) belongs to V1 , then w1(2k / 12) Å c1(k). Thus, c1(k) Å 0 for all k √ Z, and therefore w1 Å 0. Hence, we conclude that gV1 is the oblique projection of g on V1 in a direction parallel to W1 , and we denote it by gV1 Å PV1//W1g. This means that P2V1//W1 Å PV1//W1 and that the error g 0 PV1//W1g belongs to W1 . Similarly, we conclude that gW1 Å PW1//V1g is the oblique projection of g on W1 in a direction parallel to V1 . We also note that gW1 Å PW1//V1g Å g 0 PV1//W1g, and immediately get that for any g √ V0 , we have g Å PV1//W1g / PW1//V1g.

(2.4)

The algorithm for computing the coefficients (c1(k), d1(k)) of the projections gV1 Å ( c1(k)f1,k and gW1 Å ( d1(k)c1,k from the coefficients c0(k) of g Å ( c0(k)f0,k is called the analysis algorithm or decomposition algorithm (here f0,k(x) Å x[0,1)(x 0 k)). The inverse algorithm, which calculates the coefficients c0(k) from (c1(k), d1(k)), is the synthesis algorithm or reconstruction algorithm. The decomposition and reconstruc-

6111$$0211

06-10-97 11:03:28

achaa

AP: ACHA

238

AKRAM ALDROUBI

FIG. 1. Perfect reconstruction filter-bank for computing the Haar oblique wavelet transform of order zero. The transfer functions of the filters are obtained by setting z Å e0i2pf.

tion algorithms consist of a series of convolutions down-sampling and up-sampling, as depicted in Fig. 1. Using the triangle inequality and the fact that {f1,k; k √ Z} and {c1,k; k √ Z} form orthonormal bases of V1 and W1 , respectively, we deduce that \ ∑ c1(k)f1,k / k√Z

∑ d1(k)c1,k\2L2 £ 2(\c1\2l2 / \d1\2l2).

(2.5)

k√Z

Moreover, using the fact that »f1,k(x), c1,l(x)…L2 Å 201/2d0(k 0 l), we conclude that for £1 Å ( c1(k)f1,k and w1 Å ( d1(k)c1,k , we have É»£1 , w1…L2É Å 201/2É»c1 , d1…l2É £ 201/2\c1\l2\d1\l2 .

(2.6)

It follows that the angle us between any two vectors £1 √ V1 and w1 √ W1 satisfies cos(us) £ 201/2 and that \£1 / w1\2L2 § \£1\2L2 / \w1\2L2 0 2É»£1 , w1…L2É § m(\£1\2L2 / \w1\2L2),

(2.7)

where 0 õ 1 0 cos(us) Å m. The orthonormality of the bases then implies that m(\c1\2l2 / \d1\2l2) £ \ ∑ c1(k)f1,k / ∑ d1(k)c1,k\2L2 . k√Z

(2.8)

k√Z

Therefore, Identity (2.4) and Inequalities (2.5) and (2.8) imply that the set {f1,k, c1,k; k √ Z} forms a Riesz basis for the space of piecewise constants with integer knot-points V0 . For any function g √ L2 , we can choose a sufficiently small value J2 so that the orthogonal projection gJ2 Å PVJ g of g into the space of piecewise constants VJ2 is 2

6111$$0211

06-10-97 11:03:28

achaa

AP: ACHA

239

OBLIQUE MULTIWAVELETS

arbitrarily close to g. Then, we can repeatedly apply the decomposition algorithm in (2.4) to obtain J1

gJ2 Å PVJ //WJ gJ2 / ∑ PWj //VjgJ2 . 1

(2.9)

1

jÅJ2/1

It follows that, for any two integers J2 õ J1 , the set {fJ1,k , cj,k(x) Å 20j/2x[1,2)(20jx 0 2k); k √ Z, J2 / 1 £ j £ J1} is a Riesz basis of VJ2 (here fJ1,k(x) Å 20J1/2x[0,1)(20J1x 0 k)). However, by letting J2 r 0` and J1 r `, we do not get a Riesz basis of L2 . On the other hand, any finite number of functions in this set are linearly independent, and finite linear combinations of this set are dense in L2 . To summarize, we have constructed sets of functions having the following properties (i) for any j, the set {fj,k ; k √ Z} is a Riesz basis of Vj ; (ii) for any j, the set {cj,k ; k √ Z} is a Riesz basis of Wj ; (iii) for any j, Vj / Wj Å Vj01 ; (iv) for any J2 õ J1 , the set {fJ1,k , cj,k; k √ Z, J2 / 1 £ j £ J1} is a Riesz basis of VJ2 ; (v) for any J1 the set {fJ1,k , cj,k ; k √ Z, j £ J1} is dense in L2 . What is interesting is that we have a fast computational algorithm for calculating the expansion in terms of the basis {fJ1,k , cj,k ; k √ Z, J2 / 1 £ j £ J1}. Since all the convolution operators in the decomposition/reconstruction have finitely many nonzero values (see Fig. 1), the algorithm is also stable. Because of the properties (i)–(v) above, we call the set {fJ1,k , cj,k ; k √ Z, j £ J1} an inductive-limit Riesz basis of L2 , and we call the spaces Wj that are generated by the Riesz bases, {cj,k(x) Å 20j/2c(20jx 0 k); k √ Z}, oblique wavelet spaces generated by the oblique wavelet c. Unlike other wavelets, the function c0 Å 21/2x[1/2,1) does not have an average of zero. Thus, it is not an orthogonal, semiorthogonal, or biorthogonal wavelet [5] (see also [26, p. 270]). In contrast with biorthogonal wavelets, this Haar oblique wavelet has no associated dual wavelet spaces or a dual wavelet basis. However, we should note that orthogonal, semiorthogonal, or biorthogonal wavelet bases also satisfy properties (i)–(v) above and, therefore, are special cases of oblique wavelet bases. Remark 2.1. An interesting feature of the Haar oblique wavelet that we have constructed is that

6111$$0211

06-10-97 11:03:28

achaa

AP: ACHA

240

AKRAM ALDROUBI

(a) the wavelet is simply a contraction of the scaling function; (b) the support of the wavelet is half the size of the Haar orthogonal wavelet; (c) remarkably, the oblique wavelet transform is faster than the Haar transform (twice as fast in the analysis, and twice as fast in the reconstruction). Remark 2.2. It is important to note that inductive-limit Riesz bases are not necessarily bases, since the expansion of a vector is not necessarily unique. What is true, is that inductive-limit Riesz bases are complete in L2 , and that they are the union of disjoint Riesz bases. We will see that the construction of this section can be generalized for the uniwavelet and multiwavelet cases. For example, we will be able to construct supercompact polynomial spline oblique wavelet bases of order n that have supports of size n, which is half the size of the shortest wavelets that have been constructed so far (e.g., the Cohen–Daubechies–Feauveau biorthogonal wavelets of order n [10]). Of course, these new bases are not Riesz bases of L2 but inductive-limit wavelet bases of L2 , and they can be viewed as an h-type hierarchical bases (see Section 7).

3. OBLIQUE MULTIWAVELET BASES

Let {Vj}j√Z be a multiresolution of L2 generated by a multiscaling function F(x) Å (f (x), f2(x), f3(x), . . . , fr(x))T (see Definition (1.1)) and consider the subspace W1 , V0 defined by 1

H∑ ∑ r

W1 Å

iÅ1 k√Z

di1(k)201/2ci

S D

J

x 0 k ; d i1 √ l2 2

Å { ∑ D1T(k)C1(x 0 2k); D1 √ lr2},

(3.1)

k√Z

where the vector function C(x/2) Å 21/2C1(x) Å 21/2(c11(x), . . . , cr1(x))T is a linear combination of the basis {F(x 0 k), k √ Z} given by

C

SD

x Å 2 ∑ (d1∗G1)(k)F(x 0 k). 2 k√Z

(3.2)

Our goal is to find conditions that will guarantee that W1 is a wavelet space associated with our starting multiresolution. For this purpose, we consider a function g(x) Å ( C0T(k)F(x 0 k) in V0 . As we did in the Haar oblique wavelet example, we wish to decompose this function g into the sum of two functions, gV1 √ V1 and gW1 √ W1 . Equivalently, we want to find a sequence G1(k) such that C defined in (3.2) is a wavelet function associated with F. To perform this task, we consider the candidate function ga(x) Å ∑ C1T(k)F1(x 0 2k) / ∑ D1T(k)C1(x 0 2k), k√Z

6111$$0211

k√Z

06-10-97 11:03:28

achaa

AP: ACHA

(3.3)

241

OBLIQUE MULTIWAVELETS

where the coefficients C1(k) and D1(k) are obtained from the coefficients C0(k) by generalized convolutions (see definition in (1.14)) followed by a downsampling operator, C1 Å 21/2 f2[H2∗C0],

(3.4)

D1 Å 21/2 f2[d01∗G2∗C0],

(3.5)

where f2[l] is the downsampling operator which assigns to the sequence B(k) the downsampled sequence f2[B](k) Å B(2k)

∀k √ Z.

(3.6)

The mappings C0 r H2∗C0 and C0 r G2∗C0 are bounded operators on lr2 if and only ˆ 2( f )\ £ M for almost all f √ [0, ˆ 2( f )\ £ M and \G if there exists M õ ` such that \H ˆ 2( f )\ are the matrix-norms of H ˆ 2( f ), respectively ˆ 2( f )\ and \G ˆ 2( f ) and G 1], where \H (see [1, Theorem 2.2]) (we can use any matrix-norm since, for r 1 r matrices, all norms are equivalent). In this paper, the norm of lr2 is the natural norm derived from l2 ; i.e., the norm of a sequence C(k) Å (c1, . . . , cr)T in lr2 is simply given by r

\C\ 2lr2 Å

∑ \ci\2l2 . iÅ1

If we rewrite the function ga in terms of the vector function F using (3.3)–(3.5), together with relations (1.3) and (1.12), we obtain ga(x) Å 2 ∑ (F2[f2[C0T∗H2T]]∗H1)(k)F(x 0 k) k√Z

/2

∑ (F2[f2[d01∗C0T∗G2T]]∗d1∗G1)(k)F(x 0 k), (3.7) k√Z

where F2[l] is the upsampling operator which assigns to the sequence B(k) the upsampled sequence F2[B](2k) Å B(k), F2[B](2k / 1) Å 0,

∀k √ Z, ∀k √ Z.

(3.8)

Using the fact that {f1(x 0 k), . . . , fr(x 0 k); k √ Z} is a basis of V0 , the functions g and ga will be equal if and only if 2 F2[f2[C0T∗H2T]]∗H1 / 2 F2[f2[d01∗C0T∗G2T]]∗d1∗G1 Å C0T .

(3.9)

Taking the Fourier transform of (3.9), we obtain the matrix equation

[CO 0T( f ) CO 0T( f 0 12)]

6111$$0211

3

HO 2T( f )HO 1( f ) / GO 2T( f )GO 1( f ) 0 Ir

4

HO 2T( f 0 12)HO 1( f ) 0 GO 2T( f 0 12)GO 1( f )

06-10-97 11:03:28

achaa

Å 0,

AP: ACHA

(3.10)

242

AKRAM ALDROUBI

where Ir is the r 1 r identity matrix. In our computation of the Fourier transform of (3.9), we have used the fact that the Fourier transform of a generalized convolution ˆ ( f )B ˆ ( f ) and the two identities, (A∗B)(k) is the matrix product A Fourier(F2[A])( f ) Å AO (2f ),

(3.11)

S S D S DD

Fourier(f2[A])( f ) Å 201 AO

f f01 / AO 2 2

,

(3.12)

where A(k) is an arbitrary matrix-sequence whose entries [A]i,j(k) are in l2 . Equation (3.10) is satisfied for all sequences C0 √ lr2 if and only if the right-hand vector in Eq. (3.10) is identically zero: HO 2T( f )HO 1( f ) / GO 2T( f )GO 1( f ) Å Ir , HO 2T( f 0 12)HO 1( f ) 0 GO 2T( f 0 12)GO 1( f ) Å 0.

(3.13)

By taking the transpose of the two equations in (3.13), we can rewrite them in the matrix form

3

HO 1T( f ) T 1

GO 1T( f ) 1 2

1 2

T 1

4

HO ( f 0 ) 0GO ( f 0 )

F G FG HO 2( f )

Å

GO 2( f )

Ir 0

.

(3.14)

The system of functional equations above is a 2r 2 equations with the 2r 2 unknowns ˆ 2( f ). If the 2r 1 2r matrix function ˆ 2( f ) and G H

SO 1( f ) Å

F

G

HO 1T( f ) GO 1T(f) HO 1T( f 0 12) 0GO 1T( f 0 12)

(3.15)

ˆ 1( f ) are 1-periodic in f ), then ˆ 1( f ) and G is invertible (Sˆ1( f ) is 1-periodic in f since H ˆ 2( f ) of (3.14) as ˆ 2( f ) and G its inverse can be written in terms of the solutions H T SO 01 1 ( f ) Å SO 2 ( f ) Å

F

G

HO 2( f ) HO 2( f 0 12) GO 2( f ) 0GO 2(f 0 12)

.

(3.16)

Let \S\ denote the matrix-norm of a matrix S (exactly which norm is not important since, for finite dimensional spaces, all norms are equivalent). We are now ready to state ˆ THEOREM 3.1. If SO 01 1 ( f ) exists for almost all f, and if the matrix-norms of S1( f ) 01 and SO 1 ( f ) are uniformly bounded; i.e., if there exist two constants m ú 0 and M ú 0 independent of f such that \Sˆ1( f )\ £ M and \SO 01 1 ( f )\ £ m for almost all f, then the ˆ ˆ inverse Fourier transform of Hi( f ) (Gi( f )), i Å 1, 2, is an r 1 r matrix-sequence Hi(k) (Gi(k)) whose entries belong to l2 , and the mappings

6111$$0211

06-10-97 11:03:28

achaa

AP: ACHA

243

OBLIQUE MULTIWAVELETS

Hi∗l: lr2 r lr2 ,

i Å 1, 2,

(3.17)

Gi∗l: lr2 r lr2 ,

i Å 1, 2,

(3.18)

are bounded convolution operators on lr2 . Moreover, the operators R1: [C1 D1]T r 21/2 F2[C1T]∗H1 / 21/2 F2[D1T]∗d1∗G1

(3.19)

R2: [C1 D1]T r 21/2 F2[C1T]∗H2 / 21/2 F2[D1T]∗d1∗G2

(3.20)

from lr2 1 lr2 into lr2 are one–one and onto bounded linear operators; i.e., they constitute isomorphisms between lr2 1 lr2 and lr2 . The proof of this theorem, as well as those of Theorems 3.2 and 3.3 below, are presented in Section 3.1. The operators C0 r H2∗C0 and C0 r G2∗C0 (defined in Theorem 3.1) appear in the decomposition (3.3)–(3.5), while the operator defined by (3.19) gives the coefficients C0(k) in the basis of V0 from the coefficients C1(k) and D1(k) of the vector ( C1T(k)F1(x 0 2k) / D1T(k)C1(x 0 2k). Using the same assumptions as those in Theorem 3.1, we now introduce two operators, PV1//W1: V0 r V1 and PW1//V1: V0 r W1 . For a function g(x) Å ( C0T(k)F(x 0 k) in V0 , these operators are defined to be PV1//W1g(x):Å 21/2 ∑ f2[C0T∗H2T](k)F1(x 0 2k),

(3.21)

k√Z

PW1//V1g(x) :Å 21/2 ∑ f2[d01∗C0T∗G2T](k)C1(x 0 2k),

(3.22)

k√Z

where H1(k), G1(k), H2(k), and G2(k) are the matrix-sequences of Theorem 3.1, and where d01 is the impulse located at k Å 01 (see Definition (1.13)). We have the following theorem. THEOREM 3.2. The intersection V1 > W1 Å {0}. Moreover, the operators PV1//W1 and PW1//V1 are projectors (not necessarily orthogonal projectors), and for any function g √ V0 , the following decomposition holds: g Å PV1//W1g / PW1//V1g.

(3.23)

The operator PV1//W1 is the projection on the space V1 in a direction parallel to W1 , and the complementary operator PW1//V1 is the projection on the space W1 in a direction parallel to V1 . We are now ready to state THEOREM 3.3.

If the conditions of Theorem 3.1 are satisfied, then the function C1(x) Å 201/2C

6111$$0211

SD

x Å 21/2 ∑ (d1∗G1)(k)F(x 0 k) 2 k√Z

06-10-97 11:03:28

achaa

AP: ACHA

(3.24)

244

AKRAM ALDROUBI

is an oblique multiwavelet. It is associated with the multiresolution generated by the multiscaling function F, whose two-scale sequence is H1(k). The set {cij,k; k √ Z, i Å 1, . . . , r.} is a Riesz basis of Wj ; for any J1 √ Z, the set {fiJ1,k , cij,k; k √ Z, j £ J1 , i Å 1, . . . , r.} is an inductive-limit Riesz basis of L2 ; and for any J2 √ Z, the set {fiJ1,k, cij,k; k √ Z, J2 / 1 £ j £ J1 , i Å 1, . . . , r.} is a Riesz basis of VJ2 . Moreover, for any g √ VJ2 , we have the decomposition J1

g Å PVJ //WJ g / ∑ PWj//Vjg. 1

(3.25)

1

jÅJ2/1

Since V0` is dense in L2 , it follows that for any e ú 0, there exist J2 and gJ2 √ VJ2 such that \e Å g 0 gJ2\L2 £ e, and we have the decomposition J1

gJ2 Å PVJ //WJ gJ2 / ∑ PWj //VjgJ2 / e. 1

(3.26)

1

jÅJ2/1

Moreover, the algorithm that calculates the decomposition is a fast filter-bank algorithm that consists of a series of generalized convolutions, upsampling, and downsampling, as we will describe further in Section 4. 3.1. Proof of Theorems Proof of Theorem 3.1. Since \Sˆ1( f )\ £ M, it follows that there exists a constant ˆ 1( f )\ £ M1 for almost all f. Thus, the inverse ˆ 1( f )\ £ M1 and \G M1 such that \H ˆ ˆ Fourier transforms of H1( f ) and G1( f ) are matrix-sequences with entries in l2 . Using [1, Theorem 2.2] we then conclude that the convolution operators defined by C0T r C0T∗G1T and C0T r C0T∗G1T are bounded operators on lr2 . Since \SO 01 1 ( f )\ £ M, the same arguments apply to H2 and G2 . Consider the mapping R1X defined by (3.19) where X Å [C1 D1]T. If we take the Fourier transform of R1X, and use Identity (3.11), we obtain 21/2CO 1T(2f )HO 1( f ) / 21/2e0i2pfDO 1T(2f )GO 1( f ).

(3.27)

ˆ 1( f )\ £ M1 , we deduce the upper bound ˆ 1( f )\ £ M1 and \G From the upper bounds \H estimate

6111$$0211

06-10-97 11:03:28

achaa

AP: ACHA

245

OBLIQUE MULTIWAVELETS

\(R1X) O ( f )\2Rr £ Const1(\CO 1(2f )\2Rr / \DO 1(2f )\2Rr),

(3.28)

from which we deduce the inequality

* \(R X) O ( f )\ df 1

\R1X\2lr2 Å

2r R

1

0

£ Const(\C1\2lr2 / \D1\2lr2).

(3.29)

Therefore, the mapping R1 is a bounded linear operator from lr2 1 lr2 into lr2 . We now consider the problem of finding the solution to R1X Å Y, where Y is an element of lr2 . In Fourier domain, this is equivalent to solving 21/2CO 1T(2f )HO 1( f ) / 21/2e0i2pfDO 1T(2f )GO 1( f ) Å YO ( f ).

(3.30)

If we replace the variable f by f 0 12 in the equation above and if we use the fact that all the functions appearing in the equation are 1-periodic, we obtain 21/2CO 1T(2f )HO 1( f 0 12) 0 21/2e0i2pfDO 1T(2f )GO 1( f 0 12) Å YO ( f 0 12).

(3.31)

Using Definition (3.15) for Sˆ1( f ), the two equations (3.30) and (3.31) can be replaced by SO 1( f )UO ( f )XO (2f ) Å ZO ( f ),

(3.32)

ˆ ( f ) Å 201/2[Y ˆ(f) Y ˆ ( f 0 12)]T and where where Z

UO ( f ) Å

F

Ir

0

0

0i2pf

e

G

Ir

.

(3.33)

ˆ ( f ) is a unitary matrix and since by assumption Sˆ1( f ) is invertible and the Since U matrix-norm of its inverse is bounded for almost all f, we conclude that the solution ˆ (2f ) exists, and we have an upper bound of the form X \CO 1(2f )\2Rr / \DO 1(2f )\2Rr £ Const\YO ( f )\2Rr,

(3.34)

from which it follows that (\C1\2lr2 / \D1\2lr2) £ Const\Y\2lr2 . exists and is a bounded linear operator from lr2 into lr2 1 lr2 . If we Therefore, R01 1 interchange the index 1 by the index 2 in the proof, we also deduce that R01 2 exists and is a bounded linear operator from lr2 into lr2 1 lr2 . j

6111$$0211

06-10-97 11:03:28

achaa

AP: ACHA

246

AKRAM ALDROUBI

Proof of Theorem 3.2. To show that PV1//W1 is a projector on V1 , we consider the action of PV1//W1 on a function £1(x) Å

∑ C1T(k)F1(x 0 2k)

(3.35)

k√Z

that belongs to V 1 . Using the two-scale relation (1.3) and the fact that F1(x) Å 201/2F(x/2) (see (1.2)), we write £1 in terms of the basis of V0 and obtain £1(x) Å 21/2

∑ (F2[C1T]∗H1)(k)F(x 0 k), k√Z

where F2[l] is the upsampling operator defined in (3.8). Using the equation above and the definition of PV1//W1 given by (3.21), we get that PV1//W1£1(x) Å 2 ∑ F2[C1T∗f2[H1∗H2T]](k)F1(x 0 k),

(3.36)

k√Z

where f2[l] is the downsampling operator defined by (3.6). The upsampling operator appears in the above equation because we are shifting F1 on the integers, instead of on the even integers, as in Definition (3.21) of PV1//W1 . Taking the Fourier transform of the above equation and using the two identities (3.11) and (3.12) yields

S

S D S DD

(PV1//W1£1) O ( f ) Å CO 1(2f ) HO 1( f )HO 2T( f ) / HO 1 f 0

1 T 1 HO 2 f 0 2 2

FO 1( f ).

(3.37)

ˆ 1( f 0 12)HO 2T( f 0 12) in Eq. (3.37) is equal to the identity The term HO 1( f )HO 2T( f ) / H matrix Ir . To see this, we simply note that since SO 1( f )SO 2T( f ) Å I2r , we also have SO 2( f )SO 1T( f ) Å I2r and SO 1T( f )SO 2( f ) Å I2r . Explicitly, this last relation can be written as the four system of r 1 r functional equations: HO 1( f )HO 2T( f ) / HO 1( f 0 12)HO 2T( f 0 12) Å Ir ,

(3.38)

GO 1( f )GO 2T( f ) / GO 1( f 0 12)GO 2T( f 0 12) Å Ir ,

(3.39)

HO 1( f )GO 2T( f ) / HO 1( f 0 12)GO 2T( f 0 12) Å 0,

(3.40)

GO 1( f )HO 2T( f ) / GO 1( f 0 12)HO 2T( f 0 12) Å 0.

(3.41)

ˆ 1( f 0 12)HO 2T( f 0 12) is simply the functional equation (3.38), Thus, HO 1( f )HO 2T( f ) / H and we conclude that Eq. (3.37) reduces to (PV1//W1£1) O ( f ) Å CO 1(2f )FO 1( f ).

6111$$0211

06-10-97 11:03:28

achaa

(3.42)

AP: ACHA

247

OBLIQUE MULTIWAVELETS

Equation (3.42) is precisely the Fourier transform of expression (3.35) for £1 . Hence, PV1//W1£1(x) Å ( C1T(k)F1(x 0 2k) Å £1(x). It follows that V1 is invariant under PV1//W1 , from which we conclude that PV1//W1 is a projector on V1 . In a similar way, we consider the action of PW1//V1 on a function w1(x) Å ∑ D1T(k)C1(x 0 2k) k√Z

which belongs to W1 . Using the relation (1.12) and the fact that C1(x) Å 201/2C(x/2), we write w1 in terms of the basis of V0 and obtain w1(x) Å 21/2 ∑ (F2[D1T]∗d1∗G1)(k)F(x 0 k).

(3.43)

k√Z

Using the equation above and the definition of PW1//V1 given by (3.22), we conclude that PW1//V1w1 (x) Å 2 ∑ F2[D1T∗f2[G1∗G2T]](k)C1(x 0 k).

(3.44)

k√Z

Taking the Fourier transform of the equation above and using Identity (3.39), we obtain (PW1//V1w1) O (f) Å D1(2f )CO 1( f ), from which we deduce that PW1//V1w1(x) Å ( D1T(k)C1(x 0 2k) Å w1(x). Thus, W1 is invariant under PW1//V1 , and PW1//V1 is a projector on W1 . Finally, assume that e1 √ V1 > W1 . Viewing e1 as a vector in V1 , we can write e1(x) Å ∑ E1T(k)F1(x 0 2k) k√Z

Å 21/2

∑ (F2[E1T]∗H1)(k)F(x 0 k).

(3.45)

k√Z

Using the definition of PW1//V1 , we get that PW1//V1e1(x) Å 2 ∑ F2[E1T∗f2[H1∗G2T]](k)C1(x 0 k). k√Z

Taking the Fourier transform of the equation above and using Identity (3.40), we deduce that (PW1//V1e1) O ( f ) Å 0. Thus, PW1//V1e1 Å 0. Since e1 is also in W1 and since PW1//V1 is a projector on W1 , we have that PW1//V1e1 Å e1 Å 0. Therefore, V1 > W1 Å {0}. j

6111$$0211

06-10-97 11:03:28

achaa

AP: ACHA

248

AKRAM ALDROUBI

Proof of Theorem 3.3.

Let £1 and w1 be the two functions

∑ C1T(k)F1(x 0 2k),

£1(x) Å

(3.46)

k√Z

w1(x) Å ∑ D1T(k)C1(x 0 2k),

(3.47)

k√Z

belonging to V1 and W1 , respectively, and consider the sum g Å £1 / w1 , which necessarily belongs to V0 . The L2-norm of g is given explicitly by \g\2L2 Å \£1\2L2 / \w1\2L2 / »£1 , w1…L2 / »w1 , £1…L2 .

(3.48)

We need to compute each term of the right-hand side of (3.48). Taking the Fourier transform of (3.46) and using Parseval’s identity, we obtain \£1\2L2 Å

* CO (2f )FO ( f )FO *( f )(CO (2f ))*df, T 1

1

T 1

1

(3.49)

R

where, as before, ‘‘∗’’ denotes the Hermitian transpose operator. If we perform a simple change of variable f * Å 2f, split the integral over intervals of length one, use ˆ 1( f ), the two-scale relation (1.3), and the function A ˆ ( f ) defined the 1-periodicity of C in (1.7), we then obtain

* CO ( f )AO ( f )(CO ( f ))*df, 1

\£1\2L2 Å

T 1

T 1

1

(3.50)

0

ˆ 1( f ) is given by where the matrix function A

S S D S D S D S D S D S DD

AO 1( f ) Å HO 1

f f f f01 f01 f01 AO HO * / HO 1 AO HO * 1 1 2 2 2 2 2 2

.

(3.51)

ˆ ( f ); i.e., A ˆ 1( f ) Å A ˆ ( f ). This fact follows ˆ 1( f ) is another expression of A In fact, A from Poisson’s formula, relation (3.6), and the fact that »fi1(x), fj1(x 0 2k)…L2 Å »fi(x), fj(x 0 k)…L2 .

The calculations for the terms \w1\2L2 , »£1 , w1…L2 , and »w1 , £1…L2 in (3.48) yield

* YO ( f )FO ( f )YO *( f )df, 1

\g\2L2 Å

1

1

(3.52)

1

0

ˆ 1( f ) Å [CO 1T( f ) DO 1T( f )] and where Fˆ1( f ) is the product of matrix-functions where Y given by

6111$$0211

06-10-97 11:03:28

achaa

AP: ACHA

249

OBLIQUE MULTIWAVELETS

S D S D S DS S DD S S DD

FO 1( f ) Å UO

0f T f f SO 1 QO 2 2 2

*

f 2

SO 1T

0f 2

UO

*

,

(3.53)

ˆ ( f ) is given by ˆ ( f ) is the unitary matrix defined by (3.33) and Q in which U

QO ( f ) Å

F

G

AO ( f )

0

0

AO ( f 0 12)

.

(3.54)

By assumption, F is a scaling function; therefore, we know from (1.7)–(1.8) that ˆ ( f ) is self-adjoint and strictly positive. ˆ ( f ) is strictly positive. This fact implies that Q A ˆ ( f ) is unitary and since \Sˆ1( f )\ £ ˆ Therefore, F1( f ) is self-adjoint. Moreover, since U 01 M and \SO 1 ( f )\ £ m, we conclude that there exist two constants Const1 ú 0 and Const2 ú 0 such that

* YO ( f )YO *( f )df £ * YO ( f )FO ( f )YO *( f )df £ Const * YO ( f )YO *( f )df. 1

Const1

1

1

1

0

1

1

1

0

1

2

1

1

(3.55)

0

By taking the inverse Fourier transform of (3.55), we deduce that m1(\C1\2lr2 / \D1\2lr2) £ \£1 / w1\2L2 £ M1(\C1\2lr2 / \D1\2lr2).

(3.56)

Obviously, the inequalities above hold if we set £1 Å 0, or equivalently, if we set C1 Å 0. Therefore, the set {c11(x 0 k), . . . , cr1(x 0 k)}k√Z is a Riesz basis of W1 . By combining this fact with Theorem 3.2, we see that C is an oblique multiwavelet. Since the spaces Wj are simply scaled version of W1 , it follows that {cij,k ; k √ Z, i Å 1, . . . , r.} is a Riesz basis of Wj . This property and the facts that V` Å {0}, that V0` is dense in L2 , and that Vj / Wj Å Vj01 (as established in Theorem 3.2) imply that the set {fiJ1,k , cij,k ; k √ Z, j £ J1 , i Å 1, . . . , r.} is an inductive-limit Riesz basis of L2 . The fact that {fiJ1,k, cij,k; k √ Z2, J2 / 1 £ j £ J1 , i Å 1, . . . , r.} is a Riesz basis of VJ2 is simply the repeated application of the arguments that led to Inequality (3.56), starting from a function g √ VJ2 . Finally, Identity (3.25) follows from the repeated application of (3.23) in Theorem 3.2. j 3.2. Biorthogonal Multiwavelets Consider the sequence HH 1(k) Å H2Û(k) Å H2(0k),

6111$$0211

06-10-97 11:03:28

achaa

(3.57)

AP: ACHA

250

AKRAM ALDROUBI

˜ 1(k) generates a multiscaling function FH . The sampled cross-correlaand assume that H tion matrix-sequence Xf,fH (k) between F and FH can then be defined by its entries [Xf,fH ]i,j(k) as [Xf,fH ]i,j(k) Å »fi(x), fH j(x 0 k)…L2 .

(3.58)

Using the two-scale relations for F and FH and the fact that »fi1(x), fH j1(x 0 2k)…L2 Å »fi(x), fH j(x 0 k)…L2 ,

we find the relation Xf,fH (k) Å 2f2[H1∗Xf,fH ∗H2T](k) Å 2f2[H1∗Xf,fH ∗(HH 1Û)T](k),

(3.59)

˜ 1(0k). An important feature of the identity above is that where, as before, HH 1Û(k) Å H if the cross-correlation Xf,fH has the property that Xf,fH (k) Å d0(k)Ir , then relation (3.59) is simply the inverse Fourier transform of Eq. (3.38). We now consider the bounded linear operator X r 2f2[H1∗X∗H2T]

(3.60)

that acts on r 1 r matrix-sequences with entries in l2 . In this context, Relation (3.59) simply states that the operator (3.60) has an eigenvalue l Å 1 and an associated eigenvector given by Xf,fH . Moreover, the inverse Fourier transform of (3.38) states that the impulse matrix-sequence d0(k)Ir is an eigenvector for the eigenvalue l Å 1. Therefore, if the operator (3.60) has an eigenspace of dimensionality 1 for the eigenvalue l Å 1, then the cross-correlation is, in fact, the impulse-sequence Xf,fH (k) Å d0(k)Ir . This is equivalent to stating that the multiscaling functions F and FH generating the two multiresolutions Vj and V˜j are biorthogonal, i.e., the bases of V0 and V˜0 satisfy »fi(x), fH j(x 0 k)…L2 Å d0(k)d0(i 0 j).

Moreover, if we choose the wavelet CH 1 to be CH 1(x) Å 201/2CH

SD

x Å 21/2 ∑ (d1∗GH 1)(k)FH (x 0 k), 2 k√Z

˜ 1(k) Å G2Û(k) Å G2(0k), then the cross-correlation Xc,cH (k) between the two where G multiwavelets C and CH is given by Xc,cH (k) Å 2f2[G1∗Xf,fH ∗(GH 1Û)T](k) Å 2f2[G1∗G2T](k).

6111$$0211

06-10-97 11:03:28

(3.61)

achaa

AP: ACHA

251

OBLIQUE MULTIWAVELETS

FIG. 2. Perfect reconstruction filter-bank for computing the oblique multiwavelet transform. The coefficients C0 are decomposed into two sequences C1 and D1 (left filter-bank pair). The two sequences C1 and D1 can then be combined by the right pair of filters to reconstruct C0 .

The right-hand side of the last equality is the inverse Fourier transform of (3.39). Therefore, Xc,cH (k) Å d0(k)Ir , and we conclude that C and CH are biorthogonal. ˜ 1(k) generates a multiscaling function and Therefore, under the assumptions that H that the eigenspace dimensionality is one, we can use the result in [1, Section 3] to obtain the biorthogonal decomposition of a function g √ V0 as g Å PV1⊥VH 1g / PW1⊥WH 1g,

(3.62)

where PV1⊥VH 1 is the projection on V1 in a direction orthogonal to V˜ 1 and where ˜ 1 . The algorithm for PW1⊥WH 1 is the projection on W1 in a direction orthogonal to W computing the projections is a fast filter-bank algorithm depicted in Fig. 2. In particular, if g(x) Å ( CT0(k)F(x 0 k), PV1⊥VH 1g(x) Å ( CT1(k)F1(x 0 2k), and PW1⊥WH 1g (x) Å ( D1T(k)C1(x 0 2k), then we can compute the coefficients C1(k) and D1(k) by the fast algorithms C1 Å 21/2 f2[H2∗C0], D1 Å 21/2 f2[G2∗C0]. Conversely, the coefficients C0(k) can be computed from C1(k) and D1(k) by the fast algorithm C0T Å 21/2 F2[C1T]∗H1 / 21/2 F2[D1T]∗d1∗G1 .

6111$$0211

06-10-97 11:03:28

achaa

(3.63)

AP: ACHA

252

AKRAM ALDROUBI

A function g √ L2 , can be approximated at any level J2 by the projection gJ2 Å PVJ ⊥VH J g. 2

2

Therefore, we have the decomposition J1

gJ2 Å PVJ ⊥VH J g / ∑ PWj⊥WH jgJ2 1

1

jÅJ2/1

`

Å

∑ PWj⊥WH jgJ2 .

(3.64)

jÅJ2/1

By letting J2 tend to infinity, we can recover g √ L2 by the formula `

g Å ∑ PWj⊥WH jg.

(3.65)

jÅ0`

From (3.65) it follows that for the biorthogonal case, the set {cij,k; (j, k) √ Z2, i Å 1, . . . , r.} is not only an inductive-limit Riesz basis, but also a Riesz basis of L2 . 4. PERFECT RECONSTRUCTION AND PERFECT DEMULTIPLEXING FILTER BANKS

4.1. Perfect Reconstruction Filter Banks A perfect reconstruction filter bank (PRFB) [34, 35] is depicted in Fig. 2. It consists of four operators on lr2 that provide useful, simple, and efficient means to decompose a sequence C0 in lr2 into two sequences C1 and D1 in lr2 1 lr2 , and to recover the original signal C0 from C1 and D1 . Specifically, a perfect reconstruction filter bank consists of: • An analysis component (see Fig. 2), which consists of two convolution and upsampling operators that map a sequence C0 √ lr2 into the two sequences (C1 , D1) √ lr2 1 lr2

C1 Å 21/2 f2[H2∗C0],

(4.1)

D1 Å 21/2 f2[G2∗C0].

(4.2)

• A synthesis component, which maps a set of two sequences (C1 , D1) √ lr2 1 lr2 into a sequence C0 √ lr2 by the formula

C0T Å 21/2 F2[C1T]∗H1 / 21/2 F2[D1T]∗d1∗G1 .

(4.3)

The analysis and synthesis components are not independent, since the synthesis operator is the inverse of the analysis operator. This dependency imposes some relationships between H1 , G1 , H2 , and G2 . Specifically, in the Fourier domain, they must satisfy

6111$$0211

06-10-97 11:03:28

achaa

AP: ACHA

253

OBLIQUE MULTIWAVELETS

FIG. 3. Perfect demultiplexing filter bank. Two vector-sequences Cj/1(k) Å (c1j/1(k), . . . , crj/1(k))T and Dj/1(k) Å (d1j/1(k), . . . , drj/1(k))T are mixed together to form a single vector-sequence (Cj(k) Å c1j (k), . . . , crj(k))T (left filter-bank pair). The two sequences Cj/1 and Dj/1 can then be resolved from Cj by the pair of filters on the right.

the system of Eq. (3.13). Therefore, the decomposition of a function g in terms of an oblique multiwavelet basis is the repeated application of a PRFB algorithm. In signal/image processing applications, the starting point is not the function g Å T C0 ∗F, but instead a discrete sequence s(k) which is obtained from the samples of some analog function s(t). The first step is then to find the coefficients C0(k) that correspond to the sequence s(k). One way to do this, is to assume that s(k) has been obtained by sampling the function s(t) Å (C0T∗F)(t) at the sampling points t Å k/r (r is the number of scaling functions). This assumption allows us, in general, to find the coefficients C0(k) by a simple filtering of the data s(k) [37]. Another solution is to use a prefiltering operation that corresponds to the orthogonal projection of s(t) (or a model of s(t)) onto the space V0 [36]. The attractive features of this last approach are that a solution always exists (unlike the interpolation approach described in [37]) and that we obtain the exact function if it is already contained in V0 (as is assumed in [37]). Moreover, this last approach is more realistic since the samples s(k) that we obtain from s(t) are usually not related to the chosen multiwavelet transform. Instead, these samples are related to the impulse response of the acquisition device, which can provide a model for the signal (see [36]).

4.2. Perfect Demultiplexing Filter Banks We define a perfect demultiplexing filter bank (or perfect resolving filter bank) as the set of four operators depicted in Fig. 3. This is the dual of the perfect reconstruction filter bank, and it consists of:

6111$$0211

06-10-97 11:03:28

achaa

AP: ACHA

254

AKRAM ALDROUBI

• A multiplexing operator that maps and intertwines a pair of sequences (C1 , D1) √ lr2 1 lr2 into a single sequence C0 in lr2 , as given by Eq. (4.3). • A demultiplexing operator that uses Eq. (4.1) to transform a sequence C0 √ lr2 back into the two sequences (C1 , D1) √ lr2 1 lr2 , as given by the equations in (4.1).

Since we want the multiplexing and demultiplexing operators to be inverses of each other, a simple calculation shows that H1 , G1 , H2 , and G2 must be related by the Fourier relation (3.38)–(3.41). But under our assumptions (i.e., the conditions of Theorem 3.1) these are equivalent to the perfect reconstruction conditions in (3.13), as shown in the proof of Theorem 3.2. Therefore, under the conditions of Theorem 3.1, a filter bank is a perfect reconstruction filter bank if and only if its dual is a perfect demultiplexing filter bank. When r Å 1, the equivalence between the four equations in (3.38)–(3.41) and the equations in (3.13) is well known [30, 33], although no interpretation had been given to this equivalence in the past. In contrast to our present approach, previous constructions of orthogonal and biorthogonal wavelets started by imposing one or more of the equations in (3.38)–(3.41). For example, Mallat starts by imposing the QMF condition, which is Eq. (3.38), with ˆ 2( f ) Å HO 2( f ) (HO 2( f) is the complex conjugate of H ˆ 2( f )). For the biorthogor Å 1 and H nal construction of Cohen, Daubechies, and Feauveau, the construction also starts by imposing Eq. (3.38) with r Å 1. It should also be noted that for r Å 1, the present method is dual to the lifting scheme [31a]. 5. CHANGING BASES

Let L(k) be an r 1 r matrix-sequence whose entries [L]i,j(k) are sequences in l2 . If the following conditions are satisfied \LO ( f )\ £ M

a.e.,

(5.1)

ˆ 01( f ) exist for almost all f and L \LO 01( f )\ £ m,

(5.2)

then, as a consequence of [1, Theorem 2.2], we have that FÉ(x) Å

∑ L(k)F(x 0 k)

is a multiscaling function for the same spaces Vj generated by F. Sequences L(k) satisfying the properties above will be called admissible. A simple calculation shows that the generating sequence H1É(k) associated with FÉ is related to the generating sequence H1(k) associated with F by the relation H1É Å F2[L]∗H1∗L01.

(5.3)

A similar result holds for wavelets. Thus, if M(k) is an r 1 r matrix-sequence satisfying ˆ ( f )\ £ M and M ˆ 01( f ) exist for almost all f and \M ˆ 01( f )\ £ m, then CÉ(x) Å \M

6111$$0211

06-10-97 11:03:28

achaa

AP: ACHA

255

OBLIQUE MULTIWAVELETS

( M(k)C(x 0 k) is an oblique multiwavelet function generating the same spaces Wj as C. The sequence G1É associated with CÉ is related to the sequence G1 associated with C by the relation

G1É Å F2[M]∗G1 .

(5.4)

5.1. Orthogonal Bases Let A(k) be the sampled autocorrelation function of F, as defined by (1.6). Since ˆ ( f ) is positive and we are assuming that F is a multiscaling function, we know that A o 01/2 ˆ ˆ self-adjoint for almost all f. Thus, we define L ( f ) Å A ( f ). With this choice for Lo, we have that Lo(k) is a well-defined admissible matrix-sequence. Moreover, a simple calculation shows that the sampled autocorrelation of the new multiscaling function Fo(x) Å ( Lo(k)F(x 0 k) is such that Ao(k) Å d(k)Ir . Thus, the set {fo1(x 0 k), . . . , for(x 0 k)}k√Z is an orthonormal basis of V0 , and the or sets {fo1 j,k (x), . . . , fj,k(x)}k√Z are orthonormal bases of Vj . Note that many other shiftinvariant orthogonal bases of Vj exist. 5.2. Equivalent Filter Banks We say that Sˆ1( f ) (see Definition 3.15) is equivalent to SO 1É( f ) if and only if Sˆ1( f ) and SO 1É( f ) generate identical multiresolution and wavelet spaces {Vj ; j √ Z} and {Wj ; j √ Z}. Let L(k) and M(k) be two matrix-sequences satisfying the conditions of admissibility (5.1) and (5.2), and consider the matrix SO 1É( f ) Å TO ( f )SO 1( f )RO (2f ),

(5.5)

ˆ ( f ) and T ˆ ( f ) are defined as where R RO T( f ) Å

F

G

LO ( f )

0

0

MO ( f )

(5.6)

and TO T( f ) Å

F

LO 01( f ) 0

0

G

LO ( f 0 12) 01

.

(5.7)

ˆ ( f ) and T ˆ ( f ) are invertible. Clearly, from our assumptions on L and M, the matrices R Moreover, their matrix-norms, as well as the matrix-norms of their inverses, are bounded. Therefore, the matrix SO 1É( f ) satisfies the conditions of Theorem 3.1. In fact, from the discussion of the previous section, SO 1É( f ) corresponds to the same spaces Vj

6111$$0211

06-10-97 11:03:28

achaa

AP: ACHA

256

AKRAM ALDROUBI

and Wj generated from SO 1( f ), but the bases are now generated by the new multiscaling and wavelet functions FÉ(x) Å

∑ L(k)F(x 0 k),

CÉ(x) Å

∑ M(k)F(x 0 k).

Therefore we have THEOREM 5.1. The two matrix-sequences S1(k) and S1É(k) are equivalent if and only if there exist two admissible matrix-sequences L(k) and M(k) such that Identity (5.5) is satisfied. Clearly, changing bases does not change the projection operators. However, the filter-bank algorithms differ for different bases, and a particular basis may be more useful than others for a given application. 6. CONSTRUCTION OF OBLIQUE MULTIWAVELET BASES

6.1. The Haar Oblique Wavelets We start by reviewing our example in Section 2 using the theory developed in Section 3. For this case, r Å 1, and H1(k) Å 201(d0(k) / d1(k)). For G1(k) we first choose G1(k) Å 201/2d0(k). The 2 1 2 matrix function Sˆ1( f ) is given by

SO 1( f ) Å

1 / z01 2

201/2 ,

1 0 z01 0201/2 2

(6.1)

where z Å ei2pf. We have that det(Sˆ1( f )) Å 0201/2 for all f ; thus, Sˆ1( f ) has an inverse. Moreover, it is easy to see that the conditions of Theorem 3.1 on the norms of Sˆ1( f ) and its inverse are satisfied. Therefore, Eq. (3.24) gives us the oblique wavelet c(x/2) Å x[1,2)(x). The filter bank implementing the decomposition and reconstruction algorithms is depicted in Fig. 1. If, instead of the previous choice for G1 , we set G1(k) Å 201(d01(k) 0 d0(k)), we obtain the well-known orthogonal Haar wavelet. It is easy to check that, in this case, Sˆ1( f ) is given by SO 1( f ) Å 201

F

G

1 / z01 z 0 1

(6.2)

1 0 z01 z / 1

and that det(Sˆ1( f )) Å 1 for all f. As expected, the wavelet function given by (3.24) is precisely the orthogonal Haar wavelet. ˆ 1( f ) ˆ 1( f ) is to take the trigonometric polynomial given by G Another choice for G n Å (z 0 1) . For n Å 0, 1, we get the previous two cases. For n Å 2 we get that

6111$$0211

06-10-97 11:03:28

achaa

AP: ACHA

257

OBLIQUE MULTIWAVELETS

det(Sˆ1( f )) Å 03 0 e0i2pf, which is nonzero for all f. Hence, we again obtain oblique wavelets, but now with two vanishing moments. 6.2. Multiwavelets Based on the Cohen–Daubechies–Plonka Multiscaling Function The Cohen – Daubechies – Plonka multiscaling function described in [9] has the 2 1 2 generating sequence, which can be defined in terms of z Å e0i2pf by HO 1( f ) Å

F

G

1 2 3 1 3 5 1 z / 2z / z (z 0 2z / z )/2 . 4 1/32 z

(6.3)

We select sequences Gn1 of the form GO n1( f ) Å (z 0 1)nI2 . For n Å 0, we easily check that the determinant of SO 01( f ) is zero for all f. Therefore, SO 01( f ) is not invertible and we cannot find a multiwavelet for this choice (the superscript n in Sn1 refers to the choice n for Gn1). For n Å 1, the determinant of SO 11(z) is given by det(SO 11(z)) Å

z2 (63z4 / 194z2 0 1), 256

(6.4)

which is nonzero for z Å e0i2pf. We conclude that the inverse of SO 12( f ) exists, and its entries are ratios of trigonometric polynomials. It follows that the conditions of Theorem 3.1 are satisfied. Therefore, the choice GO 11( f ) Å (z 0 1)I2 gives rise to an oblique multiwavelet. It should be noted that this multiwavelet has its moment of order 1 equal to zero. Using Mathematica, we have checked that n Å 3, 5, 7 also gives rise to multiwavelets. For n Å 2 and, in fact for any even integer n Å 2m, a simple calculation shows 1 2m that SO 2m 1 ({4 ) is singular. Since SO 1 ( f ) is continuous in f, the matrix-norm of the inverse 2m of SO 1 ( f ) cannot be bounded for almost all f. Thus, it is not possible to construct 2m oblique multiwavelets using sequences of the form GO 2m 1 ( f ) Å (z 0 1) I2 . 6.3. Some Simple Choices of Oblique Multiwavelets A first candidate for an oblique multiwavelet is to choose G1(k) Å d0(k)Ir , or ˆ 1( f ) Å Ir . If the conditions of Theorem 3.1 are satisfied, then we equivalently, G obtain an oblique multiwavelet given by Eq. (3.24). This is the simplest multiwavelet. It has the simplest synthesis filter implementation. Moreover, if the functions forming the multiscaling function F have compact support (equivalently, if H1(k) has compact support), then C will have compact support, as well. In fact, the support of C will be half the size of the support of F. The analysis filters H2 and G2 will not have compact support, in general. However, the Fourier transforms of their entries are ratios of trigonometric polynomials. Therefore, it is possible to implement them exactly, using recursive filtering algorithms similar to the ones described in [32]. From this perspective, these filters can be considered as if they were compactly supported. The type of oblique multiwavelets generated by this method will not generate Riesz bases of L2 in general; instead, they generate Riesz bases for Wj and inductive-limit bases of L2 that are associated with a fast filter-bank algorithm. To obtain Riesz bases

6111$$0211

06-10-97 11:03:28

achaa

AP: ACHA

258

AKRAM ALDROUBI

of L2 , the possible candidates should have zero moments. For this purpose, G1 can ˆ 1( f ) Å (1 0 ei2pf)nIr . be chosen to be of the form G 1 ˆ 1( f ) Å HO * Another interesting candidate is given by G 1 ( f 0 2). For r Å 1 and for the case where H1 satisfies a QMF condition, this choice results in an orthogonal wavelet basis. However, r x 1 will not give rise to an orthogonal basis, in general. 6.4. General Construction Since we are dealing with real scaling and wavelet functions, our two-scale generating sequences H1(k) are real sequences, and we wish to construct sequences G1(k) that are real, as well. For this reason, we must also require that in Fourier domain, ˆ 1(0f ) Å GO 1( f ). This requirement is always satisfied for the two simple constructions G in the previous section. ˆ 1( f0) and G ˆ 1( f0 0 12) such In general, for each f0 in the interval [0, 14], we choose G ˆ 1( f )H ˆ 1( f 0 12)] must have rank r, and in this case, there that Sˆ1( f ) has an inverse ([H ˆ 1( f ) in the intervals are infinitely many possible choices). This fixes the values of G 1 1 1 1 1 [02, 04] < [0, 4]. The remaining part of the interval [02, 2] can be computed by using ˆ 1( f ) Å GO 1(0f ). Because H1(k) is real, we have that H ˆ 1( f ) Å HO 1(0f ). the formula G ˆ 1( f ) is 1-periodic, we conclude that our construction Therefore, using the fact that H is such that Sˆ1( f ) Å SO 1(0f ), that SO 01 1 ( f ) exists for almost all f and that the oblique wavelet corresponding to our choice of G1 is real.

7. OBLIQUE WAVELETS AND HIERARCHICAL BASES

7.1. The Piecewise Linear Hierarchical Basis For the case of the piecewise linear function spaces Vj , the generating sequence H1(k) Å 401(d0(k) / 2d1(k) / d2(k)) gives rise to the scaling function w(x) Å b1(x), which is the B-spline (or hat) function of order one, if we choose G1(k) Å d1(k), then the 2 1 2 matrix function Sˆ1( f ) is given by

SO 1( f ) Å

1 / 2z01 / z02 01 z 4 1 0 z01 / z02 4

, z

(7.1)

01

where z Å ei2pf. We have that det(Sˆ1( f )) Å e0i4pf; thus, Sˆ1( f ) has an inverse. Moreover, it is easy to see that the conditions of Theorem 3.1 on the norms of Sˆ1( f ) and its inverse are satisfied (for this case, m Å M Å 1). Therefore, Eq. (3.24) gives us the oblique wavelet c(x/2) Å b1(x 0 2), which, except for the shift, is the same as the scaling function that we have started with. This example is precisely the h-type hierarchical basis used in finite element [4, 38] (see Fig. 4). Unlike orthogonal, semiorthogonal, and biorthogonal wavelets, this hierarchical wavelet basis does not have an average of zero; however, it is an oblique wavelet basis.

6111$$0211

06-10-97 11:03:28

achaa

AP: ACHA

259

OBLIQUE MULTIWAVELETS

FIG. 4.

The piecewise linear hierarchical basis.

7.2. Smooth Hierarchical Basis and Supercompact Wavelets The example in the previous section suggests a simple construction of hierarchical bases. In particular, if the multiresolution spaces V 3j are the piecewise polynomial splines of order 3, generated by the B-spline scaling function b3, then we know that the Fourier transform of the generating sequence is HO 31( f ) Å 204(1 / z01)3. If we choose G1(k) Å d1(k), then the 2 1 2 matrix function Sˆ1( f ) is given by SO 1( f ) Å

F

G

204(1 / z01)3 z01 , 204(1 0 z01)3 z01

(7.2)

where z Å ei2pf. The determinant det(Sˆ1( f )) Å 203e0i4pf(3 / e0i4pf ) is nonvanishing, and all the conditions of Theorem 3.1 are satisfied. Therefore, we obtain the oblique wavelet which is a twice continuously differentiable h-type hierarchical basis having the following properties: (1) the wavelet is symmetrical, and it is a spline wavelet of order 3; (2) the size of the support of c is 2, which is shorter (by a factor 2) than the support of any spline biorthogonal (or semiorthogonal or orthogonal) wavelet of order 3. (3) the expansion of any vector in this hierarchical basis can be obtained by a perfect reconstruction filter bank. Remark 7.1. The synthesis filters H1(k) and G1(k) have finite support and their implementation has a complexity that is proportional to the size of their support. Since the analysis filters are ratios of trigonometric polynomials, the sequences H2(k) and G2(k) decay exponentially fast. Thus, for their implementation, the analysis filters can be truncated. We can also implement the analysis filters, exactly, by the iterative method developed in [32], which has the same complexity as the corresponding finitely supported synthesis filters. 7.3. The Cubic Hermite Spline Hierarchical Multiwavelet Finally, we construct a hierarchical multiwavelet by starting from the cubic Hermite spline multiscaling function which has the 2 1 2 generating sequence given by [19]

6111$$0211

06-10-97 11:03:28

achaa

AP: ACHA

260

AKRAM ALDROUBI

FIG. 5. The Hermite cubic spline oblique multiwavelet C(x) Å (c1(x), c2(x)) that generates the hierarchical Hermite cubic spline basis.

H1(0) Å

F

G

1/2 3/4 , 01/8 01/8

H1(1) Å

F G

1 0 , 0 1/2

H1(2) Å

F

G

1/2 03/4 . 1/8 01/8

(7.3)

We select sequences G1(k) to be G1(k) Å d0(k)I2 , where I2 is the 2 1 2 identity matrix. We easily check that the determinant of the 4 1 4 matrix Sˆ1( f ) is nonzero for all f. Thus, the inverse of Sˆ1( f ) exists, and its entries are ratios of trigonometric polynomials. It follows that the conditions of Theorem 3.1 are satisfied. Therefore, the oblique multiwavelet that we obtain is precisely a contraction (by a factor of 2) of the cubic Hermite spline functions that we have started with. The cubic Hermite spline hierarchical basis, which is depicted in Fig. 5, allows us to interpolate functions and their

TABLE 1 Reconstruction Filters: H1 , G1 : [H1]i, j Å hi, j , [G1]i, j Å gi, j k

h11

h12

h21

h22

g11

g12

g21

g22

01 0 1

1/4 1/2 1/4

3/8 0 03/8

01/16 0 1/16

01/16 1/4 01/16

1 01/16 0

0 0 0

0 0 0

1 1 0

6111$$0211

06-10-97 11:03:28

achaa

AP: ACHA

261

OBLIQUE MULTIWAVELETS

TABLE 2 Analysis Filters: H2 , G2 : [H2]i, j Å hi, j , [G2]i, j Å gi, j k 014 013 012 011 010 09 08 07 06 05 04 03 02 01 0 1 2 3 4 5 6 7 8 9 10 11 12

h11 0 0 0.0000047 0 0.0000468 0 0.0004634 0 0.0045871 0 0.0454077 0 0.4494897 0 0.4494897 0 0.0454077 0 0.004587 0 0.0004634 0 0.0000468 0 0.0000047 0 0

h12 0 0 00.0000039 0 00.0000382 0 00.0003786 0 00.0037454 0 00.0370752 0 00.3670068 0 0.3670068 0 0.0370752 0 0.0037454 0 0.0003786 0 0.0000382 0 0.0000039 0 0

h21

h22

0.0000023 0 0.0000232 0 0.0002293 0 0.0022701 0 0.0224721 0 0.2224513 0 2.2020410 0 02.2020410 0 00.2224513 0 00.0224721 0 00.0022701 0 00.0002293 0 00.00002317 0 00.0000023

00.0000019 0 00.0000189 0 00.0001872 0 00.0018536 0 00.0183484 0 00.1816307 0 01.797959 0 01.797959 0 00.1816307 0 00.0183484 0 00.0018536 0 00.0001872 0 00.000019 0 00.0000019

g11

g12

0 0 0 00.0000024 0 00.0000234 0 00.0002317 0 00.0022936 0 00.0227038 0 00.2247449 0.5 00.2247449 0 00.0227038 0 00.0022936 0 00.0002317 0 00.0000234 0 00.0000024 0

0 0 0 0.0000019 0 0.0000191 0 0.0001892 0 0.0018727 0

0.0185376 0 0.1835034 0 00.1835034 0 00.0185376 0 00.0018727 0 00.0001892 0 00.0000191 0 00.0000019 0

g21 0 0 0 00.0000058 0 00.0000573 0 00.0005675 0 00.0056181 0 00.0556128 0 00.5505103 0 0.5505103 0 0.0556128 0 0.0056180 0 0.0005675 0 0.0000573 0 0.0000058 0

g22 0 0 0 0.0000047 0 0.0000468 0 0.0004634 0 0.0045871 0 0.0454077 0 0.4494897 0.5 0.4494897 0 0.0454077 0 0.0045871 0 0.0004634 0 0.0000468 0 0.0000047 0

derivatives at the grid points. The filters H1 , G1 , H2 , and G2 implementing this hierarchical basis are given in Table 1 and Table 2.

ACKNOWLEDGMENTS I thank Douglas Hardin for very stimulating discussions. I am also grateful to Manos Papadakis for his careful reading of this manuscript and for his insightful comments. I also thank Richard Zalik for directing me to some very useful references, Michael Unser and Manil Suri for their helpful comments, and Emilios Dimitriadis for using his mastery to help me with Mathematica. Finally, I thank Barry Bowman for his careful editing of the manuscript.

REFERENCES 1. A. Aldroubi, Oblique projections in atomic spaces, Proc. Amer. Math. Soc. 124 (1996), 2051–2060. 2. A. Aldroubi and J. McGowan, Oblique and biorthogonal multi-wavelet bases with fast-filtering algorithms, in [22, pp. 15–26]. 3. A. Aldroubi and M. Unser, Families of multiresolution and wavelet spaces with optimal properties, Numer. Funct. Anal. Optim. 14 (1993), 417–446. 4. R. E. Bank, Hierarchical bases and the finite element method, in ‘‘Acta Numerica 1996’’ (A. Iserles, Ed.), pp. 1–43, Cambridge Univ. Press, Cambridge, 1996. 5. C. K. Chui and X. L. Shi, Bessel sequences and affine frames, Proc. Amer. Math. Soc. 1 (1993), 29– 49.

6111$$0211

06-10-97 11:03:28

achaa

AP: ACHA

262

AKRAM ALDROUBI

6. C. K. Chui and J. Z. Wang, On compactly supported spline wavelets and a duality principle, Trans. Amer. Math. Soc. 330 (1992), 903–915. 7. A. Cohen, Analyses multire´solutions et filtres miroirs en quadrature, Ann. Inst. H. Poincare´ Anal. Non Line´aire 7 (1990), 439–459. 8. A. Cohen, I. Daubechies, and J. C. Fauveau, Biorthogonal bases of compactly supported wavelets, Comm. Pure Appl. Math. 45 (1992), 485–560. 9. A. Cohen, I. Daubechies, and G. Plonka, Regularity of refinable function vectors. [preprint] 10. I. Daubechies, ‘‘Ten Lectures on Wavelets,’’ SIAM, Philadelphia, 1992. 11. I. Daubechies, Orthonormal bases of compactly supported wavelets, Comm. Pure Appl. Math. 41 (1988), 909–996. 12. C. de Boor, R. A. DeVore, and A. Ron, On the construction of multivariate pre-wavelets, Constr. Approx. 9 (1993), 123–166. 13. G. C. Donovan, J. S. Geronimo, D. P. Hardin, and P. R. Massopust, Construction of orthogonal wavelets using fractal interpolation functions. [preprint] 14. J. S. Geronimo, D. P. Hardin, and P. R. Massopust, Fractal functions and wavelet expansions based on several scaling functions, J. Approx. Theory 78 (1994), 373–401. 15. T. N. T. Goodman and S. L. Lee, Wavelets of multiplicity r, Trans. Amer. Math. Soc. 342 (1994), 307–324. 16. T. N. T. Goodman, W. S. Tang, and S. L. Lee, Wavelets wandering subspaces, Trans. Amer. Math. Soc. 338 (1993), 639–654. 17. K. Gro¨chenig, Orthogonality criteria for compactly supported scaling functions, Appl. Comput. Harmon. Anal. 1 (1994), 242–245. 18. C. Heil and D. Colella, Matrix refinement equations: existence and uniqueness. [preprint] 19. C. Heil, G. Strang, and V. Strela, Approximation order by translates of refinable functions. [preprint] 20. L. Herve`, Multi-resolution analysis of multiplicity d. Application to dyadic interpolation, Appl. Comput. Harmon. Anal. 1 (1994), 299–315. 21. B. Jawerth and W. Sweldens, An overview of wavelet based multiresolution analyses, SIAM Rev. 36 (1994), 377–412. 22. A. F. Laine and M. Unser (Eds.) ‘‘Mathematical Imaging: Wavelet Applications in Signal and Image Processing,’’ Vol. 2596, SPIE, The International Society for Optical Engineering, Bellingham, WA, 1995. 23. P. G. Lemarie´ and Y. Meyer, Ondelettes et bases hilbertiennes, Rev. Math. Iberoamer. 2 (1986). 24. S. Mallat, Multiresolution approximations and wavelet orthonormal bases of L2(R), Trans. Amer. Math. Soc. 315 (1989), 69–97. 25. J. Marasovich, ‘‘Biorthogonal Multiwavelets,’’ Dissertation, Vanderbilt University, Nashville, Tennessee, May 1996. 26. Y. Meyer, ‘‘Ope´rateurs de Calderon-Zygmund,’’ Hermann, Paris, 1990. 27. C. A. Michelli, Using the refinement equation for the construction of pre-wavelets, Numer. Algorithms 1 (1991), 75–116. 28. G. Plonka, Approximation order of shift-invariant subspaces of L2(R) generated by refinable functions. [preprint] 29. G. Plonka, Generalized spline wavelets, Constr. Approx. 12 (1996), 127–155. 30. O. Rioul, A discrete-time multiresolution theory, IEEE Trans. Signal Process. 41 (1993), 2591–2606. 31. G. Strang and V. Strela, Orthogonal multiwavelets with vanishing moments, J. Optica Eng. 33 (1994), 2104–2107. 31a.W. Sweldens, The lifting scheme: A custom-design construction of biorthogonal wavelets, Appl. Comput. Harmon. Anal. 3 (1996), 186–200. 32. M. Unser, A. Aldroubi, and M. Eden, Fast B-spline transforms for continuous image representation and interpolation, IEEE Trans. Patt. Anal. Mach. Intell. 13 (1991), 277–285.

6111$$0211

06-10-97 11:03:28

achaa

AP: ACHA

263

OBLIQUE MULTIWAVELETS

33. M. Unser, A. Aldroubi, and M. Eden, A family of polynomial spline wavelet transforms, Signal Processing 30 (1993), 141–162. 34. P. P. Vaidyanathan, ‘‘Multirate Systems and Filter Banks,’’ Prentice-Hall, Englewood Cliffs, NJ, 1993. 35. M. Vetterli and J. Kovacevic, ‘‘Wavelets and Subband Coding,’’ Prentice-Hall, Englewood Cliffs, NJ, 1995. 36. M. J. Vrhel and A. Aldroubi, Pre-filtering for the initialization of multi-wavelet transforms, in ‘‘Proc. IEEE ICASSP, April 21–24, 1997, Munich, Germany.’’ 37. X. G. Xia, J. S. Geronimo, D. P. Hardin, and B. W. Sutter, Design of prefilters for discrete multiwavelet transforms, IEEE Trans. Signal Process. 44 (1996), 25–35. 38. H. Yserentant, On the Multi-level splitting of finite element spaces, Numer. Math. 49 (1986), 379– 412.

6111$$0211

06-10-97 11:03:28

achaa

AP: ACHA