two-dimensional orthogonal wavelets with vanishing ... - CiteSeerX
TWO-DIMENSIONAL ORTHOGONAL WAVELETS WITH VANISHING MOMENTS David Stanhill and Yehoshua Y. Zeevi Department of Electrical Engineering Technion - Israel institute of technology Haifa 32000, ISRAEL fax: (972) 4 323041 email: [email protected] June 1995 CC PUB #112, (or EE PUB # 976)
Abstract We investigate a very general subset of two-dimensional, orthogonal, compactly supported wavelets. This subset includes all the wavelets with a corresponding wavelet (polyphase) matrix that can be factored as a product of factors of degree-1 in one variable. In this paper we consider in particular wavelets with vanishing moments. The number of vanishing moments that can be achieved increases with the increase of the McMillan degrees of the wavelet matrix. We design wavelets with the maximal number of vanishing moments for given McMillan degrees, by solving a set of nonlinear constraints on the free parameters de ning the wavelet matrix, and discuss their relation to regular, smooth wavelets. Design examples are given for two fundamental sampling schemes; the quincunx and the four-band separable sampling. Their relation to the, well known, one-dimensional Daubechies wavelets with vanishing moments is discussed.
This research was supported in part by the Ollendorf Center of the Department of Electrical Engineering, by the Israel Science Foundation and by the Fund for Promotion of Research at the Technion.
1
1 Introduction The M-band 2D discrete wavelet transform (DWT) is de ned by M functions, one scaling function - 0(x) and (M-1) wavelet functions - p(x), and a sampling matrix - D. For wavelets associated with multiresolution analysis the DWT can equally be speci ed by M wavelet lters - H p(!), instead of the M wavelet functions. In this work we restrict the discussion to compactly supported functions that de ne an orthogonal DWT. The orthogonal compactly supported wavelet lters are a special case of lossless FIR multirate lter banks. Since considerable work has been carried out on lossless lter banks, we adopt, in part, the notation and the terminology used in this eld. Speci cally, we limit most of the analysis to the polyphase matrix de ned by the wavelet lters (wavelet matrix) and not to the wavelet functions. Most current applications of the 2D DWT use 1D wavelet functions and lters to construct 2D wavelets, usually separable wavelets[1]. Using separable functions preserves orthogonality but imposes a severe limitation on the 2D wavelets. For example the commonly used four-band lters can not have linear phase, except for the Haar wavelet, and their orientation can not be adjusted for given speci cations. Various other methods were proposed for the design of 2D perfect reconstruction lter banks [2, 3, 4], which do not require orthogonality (paraunitarity) or nite support. The restrictions imposed on the functions/ lters still leave considerable freedom in choosing the wavelet functions/ lters. Motivated by the work of Daubchies in one dimension[5], we examine the design of 2D, orthogonal, compactly supported wavelets with maximal number of vanishing moments. In the 1D case these have become the standard wavelets, and they can be considered suboptimal in a large class of applications, such as the representation, processing and compression of 1D signals. For detailed experimental results in image compression using separable wavelets, showing the optimality of wavelets with vanishing moments, see [6]. The generalization of the 1D notion of vanishing moments to higher dimensions is straightforward, but its nature and importance is still to be determined. In 1D 2-band wavelets, the number of vanishing moments is closely related to the atness of the frequency response of the wavelet lters at the aliasing points. When dealing with more than two bands, even in 1D, the M functions are no longer uniquely de ned by each other, and the vanishing moment constraints ensure atness of the scaling lter but not of the other wavelet lters. The regularity of the wavelet functions, is also closely related to the number of vanishing moments in the 1D case[7, 8, 9]. Yet, in more than 1D no such relation has been rigorously shown, although it was suggested to hold by Daubechies and Lagarias[8]. The work done by Villemoes on the quincunx sampling scheme[10], introduces a di erent approach to the general 2D regularity measurement, and is used here to obtain lower bounds of the Holder exponent. The method used by Daubechies[5] in 1D can not be generalized to construct orthogonal wavelets with N vanishing moments in more than one dimension. The factorization of the frequency response of the scaling lter as 1 + exp(?i!) N 0 L(!); H (!) = 2 has no strait-forward generalization to more than 1D, since the fundamental theorem of 2
algebra holds only for 1D. We therefore construct 2D wavelets by rst parameterizing the space of orthogonal wavelets with compact support, and then impose a set of constraints on the parameters in order to achieve the desired number of vanishing moments. The shortcome of this method is that no complete parameterization of the space of 2D, orthogonal, compactly supported wavelets exists. In most of this work we limit the design to wavelets with wavelet matrices that are factorable as a product of 1D degree-1 factors. Such wavelets were proposed by Karlsson and Vetterli[11]. The issue of parameterization of 2D lossless FIR transfer matrices has been investigated by Venkataraman and Levy[12, 13], they have shown that almost all lossless FIR transfer matrices (wavelet matrices in our case) are factorable. Towards the completion of this manuscript our attention was drawn to a recent work by Kovacevic and Vetterli[14]. They use a construction method similar to the one proposed here, and give a few design examples. The rest of the paper is organized as follows. The rst section is devoted to the description of 2D orthogonal wavelets and lter banks in terms of their polyphase matrix. Two special subclasses are discussed, factorable and separable wavelets and lter banks. In the next section we give the conditions for obtaining wavelets with L vanishing moments, and discuss their relation to regular, smooth wavelets. Speci c design results and examples are given in the last section; rst for quincunx wavelets and then for four-band, separable sampling wavelets. A few comments on notation. The term degree refers to the polynomial degree if not otherwise speci ed. We use boldface for two-dimensional entities like x; n. We adopt the notation used in the multidimensional multirate literature; for two vectors n = (n1 ; n2)T 4 n n and z = (z1 ; z2)T we have zn = z1 z2 , and for a 2 2 matrix D we de ne, zD =4 (z1D ; z2D ; ; z1D ; z2D ; ). Capital letters are used for functions given in the z and frequency domain, such as H (z) and H (!). Finally, boldface is used also for matrix valued functions. 1
11
21
12
2
22
2 2D Wavelets and multirate lter-banks When dealing with multi-dimensional wavelets and multirate systems, the change in resolution and sampling rate is given by a matrix -D; a 2 2 matrix in the 2D case as compared with a single factor in the 1D case. For a detailed analysis of 2D sub-sampling and multirate systems see Viscito and Allebach[15], and for a more mathematical background in multidimensional wavelets see Lawton and Resniko [16]. Here we quote a few facts and de nitions relevant to the rest of this paper. For a discrete signal given on a lattice , the sub-sampled signal is de ned on the lattice D = fDn : n 2 g. The density of grid points on D is smaller by a factor of M = j det(D)j than the density on . A coset of a sub-lattice is the set of grid points on the lattice which are also on a shifted version of the sub-lattice, a coset is referred to by a vector ki , that de nes the 2D shift. There are M distinct cosets which cover the whole lattice. The points !i = 2D?T ki in the frequency plane are called points of repeated spectra or \aliasing points". We now de ne 2D wavelets. The discrete 2D wavelet transform (DWT) of a function f 2 L2 (IR2) is de ned as
Z
4 f (x) p (x) dx; (W f )pj;n = j;n
3
(1)
where
4p j p j p 2 j;n(x) = M (D x ? n) : p = 1; : : :M ? 1; j 2 ZZ ; n 2 ZZ : (M ? 1) wavelet functions - p(x) are related to the scaling function - 0(x) through
The the two-scale re nement equations:
p (x) = pM X hp (n)
(Dx ? n); (2) n where hp (n) is a 2D sequence. In this work we focus our e ort on real, compactly supported wavelets, hence the sequences hp (n) are real and of nite length. The M sequences hp (n) and their Fourier transform H p(!) can be viewed as the impulse and frequency responses of 2D FIR lters, and are called wavelet lters. H 0(!) the scaling lter is a low-pass lter and the rest are high-pass or band-pass lters. It is by means of these wavelet lters, that we obtain, the now well established, connection between wavelet theory and multirate lter banks. To assist the analysis of multirate lter banks we de ne two matrices, the polyphase and modulation matrices. In multirate systems each lter H p(!) is decomposed into M polyphase components, one for each coset
H p(!) =
MX ?1 j =0
0
e?i!kj Hjp(DT !);
(3)
where each polyphase component is given by,
Hjp(!) =4
X
hp (Dn + kj ) e?i!n : (4) n The polyphase matrix of a maximally decimated lter bank is a M M matrix, X (5) H(!) = h(n) e?i!n; n with one column for each lter and one row for each polyphase component, that is the p'th column is T p H0 (!); H1p(!); HMp ?1(!) (Note that the polyphase matrix is usually de ned as the transpose of the one de ned above, we use our de nition for future simpli cation.) The modulation matrix - HM (!) is closely related to the polyphase matrix. In this case the p'th column is built from M shifted versions of the lter H p(!), one for each aliasing frequency, that is the p`th column is T p p M H (!); H p(! + 2D?T k1 ); H p (! + 2D?T kM ?1) :
Claim 1 A maximally decimated lter bank is lossless if and only if it obeys one of the
following equivalent conditions : 1.
P hp(n)hq(Dm + n) = : p;q m; n 0
4
2. Hy(!)H(!) = I or
P hy(n + k)h(n) = I : k; n 0
3. HyM (!)HM (!) = I: Claim 2 A bank of M , 2D, FIR lters de nes a compactly supported, orthogonal 2D wavelet system if and only if : The lter bank is lossless, p Pn hp(n) = H p(! = 0) = p;0 M: The term orthogonal applies only to the wavelet lters and not to the wavelet functions, which in general are not orthogonal but constitute a tight frame with frame bounds A = B = 1. In the z-domain the polyphase matrix of FIR (causal) lters is a 2D polynomial matrix in the variables z1?1 and z2?1 . The unitary condition of the polyphase matrix is now HT (z?1)H(z) = H(z)HT (z?1) = I: (6) Polyphase matrices that obey the above condition are called paraunitary, i.e. unitary at all points jzi j = 1. We use the term wavelet matrix of rank M for the polyphase matrix of an orthogonal M -band wavelet system. Since the characteristic matrix H0 of a paraunitary polynomial matrix (PUP matrix), de ned as H0 = H(z1 = 1; z2 = 1), is a constant orthogonal matrix, any PUP matrix can always be decomposed as, H(z1; z2) = HI (z1; z2)H0; with HI (z1 = 1; z2 = 1) = I: For wavelet matrices the characteristic matrix must have a special form so as to satisfy the constraints in claim 2. These matrices are called Haar wavelet matrices[17], and are a generalization of the, well known, two-band Haar matrix. The notation - H(M ) is used for a Haar wavelet matrix of rank M . Throughout the rest of the paper we will work in the z - domain, except in section 4, where the results are better understood in the frequency domain. Most results will be stated for general PUP matrices but the design examples will be limited to wavelet matrices. The reader unfamiliar with 2D polynomial matrices, and notions such as the McMillan degrees should refer to [12], for details. We quote one result from [12] for reference and further consideration. Claim 3 A polynomial matrix is paraunitary and of McMillan degrees (m,n) if and only if; (7) det H(z1 ; z2) = z1?m z2?n :
2.1 Factorable wavelets
Denition 1 A PUP matrix of McMillan degrees (m,n) in z and z is called factorable 1
2
if it can be factored into a product of McMillan degree-1, 1D matrices, m of them in z1 and n in z2 , that is: mY +n H(z1; z2) = Hk (zk ); (8) k=1
where zk = z1 or z2 , and Hk are 1D PUP matrices of McMillan degree-1.
5
Using 1D results[18] we can write any factorable PUP M M matrix of McMillan degrees (m; n) as: mY +n (9) H(z1; z2) = [I + (zk?1 ? 1)VkVkT ]H0; k=1
where H0 is the characteristic matrix of H, Vk are unit vectors of size M , and as above, zk = z1 or z2. As mentioned earlier, for wavelet matrices H0 must satisfy an extra constraint of regularity, and is denoted by H(M ) . The number of free parameters in the above factorization (9) is ! ! M M ? 1 Nfp = + (M ? 1)(m + n); or, for wavelets + (M ? 1)(m + n); 2 2 where the binomial part comes from the characteristic (Haar wavelet) matrix and the rest from the (m + n) unit vectors of length M . In [12] it is shown that this is equal to the number of degrees of freedom of all PUP matrices of rank M and McMillan degrees (m,n), hence the set of factorable PUP matrices is not a small subset of all PUP matrices. Furthermore, in the case that one of the polynomial degrees is equal to one, all the PUP matrices are factorable[2]. In our complementary study [19] we show that a PUP matrix is factorable in the McMillan degree sense, as de ned above, if and only if it is factorable in the polynomial sense. The polynomial factorization is used for the factorization of linear phase PUP matrices.
2.2 Separable Wavelets
In most applications of the 2D DWT, separable wavelets are used. In this section it is shown how separable wavelets t into the general framework of 2D wavelets. Denition 2 We call a MN -band, 2D lter bank separable if it can be separated into a direct (Kronecker) tensor product of 1D M-band and N-band lter banks, that is (H 0(z1; z2); H 1(z1 ; z2); : : :; H MN ?1(z1 ; z2)) = (H~ 0(z1)H 0(z2); : : :; H~ 0(z1 )H N ?1(z2); H~ 1(z1 )H 0(z2); : : :; H~ M ?1(z1 )H N ?1(z2 )) = (H~ 0(z1); H~ 1(z1 ); : : :; H~ M ?1(z1 )) (H 0(z2); H 1(z2 ); : : :; H N ?1(z2)): ?1 p(z2 )gpN=0?1 are M-band and N-band 1D lter banks respectively. where fH~ p(z1)gM p=0 and fH Separable lter banks are a special case of 2D lter banks with separable sampling, i.e. with a diagonal sampling matrix. A straight forward way of de ning the polyphase matrix for separable lter banks is by choosing the coset shift vectors as: ( ! ! ! ! ! !) 0 0 0 1 1 M ? 1 ki = ; ; : : :; ; ; : : :; ; : : :; : 0 1 N ?1 0 N ?1 N?1 Hence the MN lters can be reconstructed from their polyphase components, (H 0(z); H 1(z); : : :H MN ?1(z)) = (zk ; zk ; : : : zkMN ? )H(zD ) = (1; z2; : : :z2N ?1 ; z1; : : :z1M ?1 z2N ?1 )H(z1M z2N ): 0
6
1
1
With this choice of coset shift vectors we have the following lemma. Lemma 1 A 2D MN -band PUP lter bank system is separable if and only if its polyphase matrix (with the above de ned coset shift vectors) is given by: H(z1; z2) = H~ (z1) H (z2); ~ (z1 ) and H (z2) are 1D M-band and N-band PUP matrices respectively. where H Proof: using the identity [A B ][C D] = [AC ] [BD]: (10) we have (1; z2; : : :z2N ?1 ; z1; : : :z1M ?1 z2N ?1 )H(z1M ; z2N ) = (H 0(z1; z2); : : :; H MN ?1(z1; z2)) = (H~ 0(z1); : : : H~ M (z1)) (H 0(z2 ); : : : H N (z2 )) ~ (z1M )] [(1; z2; : : :z2N ?1 )H (z2N )] = [(1; z1; : : :z1M ?1)H (z2N )] ~ (z1M ) H = [(1; z1; : : :z1M ?1) (1; z2; : : :z2N ?1 )][H ~ (z1M ) H (z2N )]; (11) = (1; z2; : : :z2N ?1 ; z1; : : :z1M ?1 z2N ?1)[H ~ (z1 )
which proves that H(z1; z2) is separable only if it can be decomposed as H(z1 ; z2) = H H(z2). The converse can be obtained by reversing the order of the equalities in (11). Using (10) once more, we get H(z1; z2) = H~ (z1) H (z2) ~ I (z1) H I (z2 )][H(0M ) H0(N )]: (12) = [H Hence, separable orthogonal MN -band lter banks have a separable characteristic matrix. We still have to show that paraunitarity is preserved. If both 1D matrices are paraunitary we have, H(z1; z2)H(z1?1; z2?1)T = H~ (z1) H (z2) H~ (z1?1) H (z2?1) T ~ ~ T ?1 T ?1 (z1 )H (z1 ) H(z2 )H (z2 ) = H = I (M ) I (N ) = I (M +N ): The converse is also true (once we x the normalization) since, A B = I (M +N ) ) A = diag(a; a; : : :; a) and B = diag(1=a; 1=a; : : :; 1=a); for some a 6= 0.
2
It is important to note that due to the fact that ~ N ? M ~ H ) = [det(H ~ I (N ))][det(I (M ) (H )] = det(H det(H ) det(H) ; from claim 3 we see that the direct tensor product of two 1D, PUP matrices of McMillan degree m and n, respectively, results in a 2D PUP matrix of McMillan degrees (Nm; Mn). 7
Theorem 1 Every 2D, MN -band, separable PUP matrix constructed from a M -band 1D PUP matrix with McMillan degree m, and a N -band 1D PUP matrix with McMillan degree n, is factorable and of McMillan degrees (Nm; Mn) Proof: Using identity (10) we have; 3 "Y # 2Y m n HI (z ; z ) = [I M + (z ? ? 1)ViViT ] 4 [I N + (z ? ? 1)Uj UjT ]5 i j 3 "Y #2 m m Y = [I M + (z ? ? 1)ViViT ] I N 4I M [I N + (z ? ? 1)Uj UjT ]5 i j n m i i Y h MN Y h MN I + (z ? ? 1)(I M Uj UjT ) = I + (z ? ? 1)(ViViT I N ) j i n m i i Y h MN Y h MN I + (z ? ? 1)Qj ; I + (z ? ? 1)Pi = 1
(
2
)
1
1
(
=1
)
2
1
=1
(
)
1
(
1
)
(
)
(
=1
)
2
1
=1
(
)
1
1
(
(
)
)
2
1
(
)
=1
=1
(
)
1
i=1
1
(
)
2
j =1
1
where Vi and Uj are unit vectors of length M and N respectively, and Pi and Qj are matrices, the nature of which is yet to be determined. De ning the vector k to be a vector of zeros except for the k'th entry which is one, we can write
Pi = ViViT I (N ) = Vi ViT
=
N X
N X k=1
k Tk
(Vi k )(Vi k )T :
k=1
(13)
Since the vectors (Vi k ) ; k = 1; : : :N are unit vectors of length MN and are orthogonal to each other, Pi = ViViT I (N ) is a rank-N symmetric projection matrix from IRMN . In a similar manner we have
Qj = I
M) U
(
T j Uj
=
M X
k=1
(k Uj )(k Uj )T ;
(14)
which again can be seen to be a rank-M symmetric projection matrix from IRMN . Summing up, we have Nm factors of McMillan degree-1 in z1 and Mn factors in z2 , thus completing the proof. It goes with out saying that the polynomial degrees of a separable PUP matrix are equal to the polynomial degree of its two 1D factors in each variable respectively. 2
3 Vanishing Moments and Regularity 1D wavelets with vanishing moments have been used extensively in the context of 1D and 2D applications of the DWT . The requirement of maximal number of vanishing moments was found to be appropriate for a variety of applications, and gave close to optimal results in 8
most cases. This success in 1D motivates us to study 2D wavelets with vanishing moments, and to apply them in image coding and representation, as well as in other 2D applications. As in the 1D case[20] we can characterize the quality of vanishing moments (or Lregularity) in several analogous manners. Theorem 2 . The following three conditions are equivalent, and can serve as a de nition of an orthogonal M-band wavelet system with L vanishing moments. 1. All moments up to order (L - 1) of the wavelet lters vanish, that is
k+q H p(!x ; !y ) X pk;q =4 hp (n)nq = (?i)k+q @ @! = 0; q xk @!y n !=0
(15)
for all p = 1; : : : (M ? 1) and k + q = 0; 1; : : :L ? 1, where nq denotes (m; n)(k;q) = mk nq . 2. All moments up to order (L - 1) of the wavelet functions vanish, that is: Z k+q ^p(!x ; !y ) @ 4 p mk;q = xq p(x) dx = (?i)k+q = 0; k q
@!x @!y
(16)
!=0
for p; k; q as above and again xq denotes xk y q .
3. The frequency response of the scaling lter has a zero of order L at the (M ? 1) aliasing frequencies, i.e. @ k+q H 0(!x ; !y ) (17) k q = 0;
@!x@!y
!=!i
for all (M ? 1) aliasing frequencies ! i and (k + q ) = 0; 1; : : :L ? 1.
Proof: 1 ) 2 : using equation (2) we can derive the recursion equation p X p Z q p mk;q = M h (n) x (Dx ? n) dx n p X p Z ? = M h (n) (D ( + n))q () dM : 0
1
n Writing (D?1 ( + n))q explicitly and using the binomial expansion we arrive at
mpk;q
! !X q q?i?j k X X k q i+j k+X 1 i q ? i = p ab i M i=0 j=0 ps;t m0(i+j?s);(k+q?i?j?t) ;
j
s=0
(18)
0
t=0
i+j s
!
!
k+q?i?j t
(19)
where D?1 = (a; b), and i = (i; j )T . Hence if pk;q = 0; 8(k + q ) < L on the right hand side of (19) we get mpk;q = 0 on the left hand side. 2 ) 1 : We use the recursion relation (19) and prove by way of induction. a) k = q = 0: p0;0 = 0 must be true for all wavelets (claim 2). 9
b) k + q = l: Given mpk;q = 0 8 (k + q ) l, we assume that pk;q = 0 8 (k + q ) < l and will prove that pk;q = 0 8 (k + q ) = l. Using the above assumption, in equation (19) we are left only with the terms in which s + t = l, that is, for k + q = l
! q k X X k = aibq?i i i j 0 Xl @X i q?i = ab
0 = mpk;q
=0 =0
i j
( + )=0
i?j
!
q p m0 j i+j;k+q?i?j 0;0 ! !1 k q A p )m00;0; i +j;l?(i+j ) i j
since m00;0 6= 0, for each p, we are left with l + 1, linear homogeneous equations for l + 1 unknowns - pi+j;l?(i+j ) . Since D is not singular the linear system is non-singular and we arrive at pi;l?i = 0 8i l. 3 ) 1 : Once again we use induction for the proof. a) k = q = 0 : p0;0 = 0 must hold for all wavelets. b) k + q = l : given that
k+q H 0(!x ; !y ) 0 = @ @! ; for all aliasing frequencies ? !i; and(k + q) l; q xk @!y !=!i
we assume that pk;q = 0; 8p = 1; : : : (M ? 1) ; k + q < l and prove that pk;q = 0 for all k + q = l. The orthogonality condition in the modulation representation (claim 1) gives
H 0(!)H p(!) +
X i
H 0(! + !i )H p(! + !i ) = p;0 ;
where the sum is over all aliasing frequencies. Di erentiating this equation l times, i.e. with (k + q ) = l, we obtain
"
@ k+q H 0(!)H p(!) + X H 0(! + ! )H p(! + ! ) i i @ !kx @ !qy i
Using condition 3,
# = 0: ! =0
! ! st q k X k?s q?t H p(!x ; !y ) X q k @ H ( ! ; ! ) @ x y 0= ; @!xs @!yt ! @!xk?s @!yq?t ! t s s t +
0
+
=0 =0
=0
=0
because of the assumption, only the term of s = t = 0 will contribute to the result, hence we are left with p p @ k+q H p(!x; !y ) M (20) = Mp = 0; p = 1; : : :; M ? 1:
!
@!xk @!yq
k;q
=0
1 ) 3: This can be shown in a very similar manner by starting from the orthogonality of the rows of the modulation matrix, that is
X j
H j (!)H j (! + !k ) = k;0 10
(21)
2
and di erentiating it.
One main di erence between the 2-band case, and M-band wavelets with M 6= 2, is that in the 2-band case the wavelet and scaling functions are fully de ned by each other. The requirement of maximum number of vanishing moments, as can be seen from the third alternative in Theorem 2, implies that the scaling lter resembles an ideal low-pass lter, and in the 2-band case, this implies that the wavelet lter is close to an ideal high-pass lter. While when we have more than one wavelet the requirement of vanishing moments does not give any information about the division of bands between the wavelet lters, and these may overlap extensively. We now deal with the issue of regularity. In the 1D case the relation between the number of vanishing moments of a wavelet function and its regularity has been investigated from several di erent viewpoints[7, 8, 9]. In all three papers cited above, it is concluded that a high number of vanishing moments is a necessary condition for obtaining high regularity, where regularity is measured by the Holder exponent. No such results exist for general N-dimensional wavelets. Furthermore, as suggested in [21, 22], for a given sampling scheme the regularity varies considerably by altering the choice of the sampling matrix and the choice of the coset vectors. The generalization of the 1D derivations is not trivial since they depend, in a crucial way, on the factorability of 1D polynomials. We face this problem, once again, in our attempts to measure the regularity, in terms of the Holder exponent, of a given wavelet function/ lter. A di erent method was suggested by Villemoes, and was adapted to the 2D quincunx case[10]. It works for Holder exponents < 1, and we use it to get a lower bound for the regularity. In the examples considered in the following section, the regularity is assessed by using a simple generalization of the technique suggested in [9]. We de ne the approximation of order-j of the scaling function 0(x) at the points x = D?j n to be : p X (22) h0j (n) = M h0 (n ? Dm)h0j?1 (m); m where h01 (n) = h0(n). The regularity of a scaling lter/function is measured by the rate of convergence of the approximation to the scaling function. Denition 3 The regularity measure of a 2D scaling function/ lter has the value , where N ? 1 < N , if it has at least N vanishing moments and ?j ?N ; 8j ; (23) max max j(l;N ?l)h0j (n)j CM lN n (
2
)
where (k;q) h0j (n) stands for the nite di erence of order (k + q ),
with and
4 (k;q)h0j (n1 ; n2) = [(k?1;q)h0j (n1 ; n2) ? (k?1;q) h0j (n1 ? 1; n2)]=M ?j=2 = [(k;q?1)h0j (n1 ; n2) ? (k;q?1) h0j (n1; n2 ? 1)]=M ?j=2
(1;0)h0j (n1 ; n2) = [h0j (n1; n2 ) ? h0j (n1 ? 1; n2)]=M ?j=2 (0;1)h0j (n1 ; n2) = [h0j (n1; n2 ) ? h0j (n1; n2 ? 1)]=M ?j=2 11
As mentioned above, it is very dicult to give general results regarding the regularity of wavelets in more than 1D. The next theorem, which is a generalization of the 1D version given in [7], and is given for a certain class of sampling matrices, shows, however, that the number of vanishing moments is related to the regularity in more than 1D as well. (Although it is stated for wavelets with compact support, it can be proven for more general wavelets with sucient decay.) Theorem 3 In the case of sampling matrices D, for which there exists an integer n such that Dn = M n I = I , a compactly supported wavelet function (x) that is orthogonal, i.e.
Z
j;n(x) i;m(x) dx = j;i n;m ;
(24)
is in C m+1 only if all its rst (m +1) moments vanish, that is, using the notation of Theorem 2 Z mk;q = xq (x) dx = 0; 8(k + q) = 0; : : :m (25)
Proof: (Based on the 1D proof [23]). First we show that the theorem is correct for m = 0. From the orthogonality, we have for every j j , Z Z 0 = (x) (Dj x ? D j ?j n ) dx = M ?j (D?j + D?j n ) () d: (26) 0
(
0)
0
0
0
If we write j = pn, and take the limit of the above expression as p grows to in nity, we have, Z Z ?pn + D?j n0 ) () d = lim ( ?p + D?j n0 ) () d 0 = plim ( D !1 p!1 =
0
Z
(D?j0 n0) () d = (D?j0 n0 )
Z
0
() d:
(27)
Note that we can change the order of the integral and the limit, since () is in C 1 and is compactly supported. This ends the proof for m = 0, since we can choose the point D?j n0 so that (D?j n0 ) 6= 0. We proceed to prove the Theorem by induction on l. Given that (x) 2 C l+1 , with l m we assume that all moments vanish for (k + q ) = i < l and prove that all the moments with (k + q ) = l, equal zero. We can use Taylor's Theorem. Since () is in C l+1 , denoting = (x; y ) we have Xl 1 Xi i ! @ i (D?j n0) ?pn (k;i?k) ? pn ? j (D ) + Rl (D?pn + D?j n0 ); (D + D n0) = i! k i ? k @x @y i=0 k=0 k (28) where the remanent term - Rl is uniformly bounded and ?j lim Rl ( + D n0) = 0 : 0
0
0
0
0
0
!0
jD?pnjl
Hence, equation (26) can be written as ! i ?j Z i Xl 1 X i @ (D n0 ) (D?pn )(k;i?k) () d 0 = @xk @y i?k i=0 i! k=0 k 0
+
Z
Rl(D?pn + D?j n0 ) () d: 0
12
(29)
Using Dn = I , we have
(D?pn )(k;i?k) = ?pi xk y i?k ; and D?pn i = ?pi jji :
According to our assumption, only the term with i = l is non zero. Hence, by dividing both terms of (29) by ?pl , we have l X l 0 = l1! k=0 k
!
@ l (D?j n0 ) Z xk yl?k () d + Z Rl ( ?pi + D?j n0 ) jjl () d: (30) @xk @y l?k j ?pi jl 0
0
Once more we take the limit of the above equation as p ! 1, and change the order of the integration and limit, to obtain
!
l X l @ l (D?j n0) Z 0 = l1! @xk @y l?k k=0 k ! l ?j Z l X l 1 @ (D n0) = l! k=0 k @xk @y l?k
Z
0
xk yl?k () d +
0
xk yl?k () d + 0:
Rl( lim p!1
?pi + D?j0 n0 ) jjl () d l ? p j j
(31)
In order to complete the proof we still have to show that each term in the above sum is zero. Since (31) is true for every point D?j n0, i.e. for every point in a set dense in IR 2, l this can be done by showing that the l + 1 functions @x@ k @y(x;yl?)k are independent. if 0
0
Xl @ l (x; y) a ;
k=0
k @xk @y l?k
we can Fourier transform this equation and integrate by parts, using the fact that all functions vanish at in nity, we have 0 0
Xl Z 1 @ l (x; y) ?i x!x e a
k=0 Xl
k=0
k
?1 @xk @y l?k
ak !xk !yl?k
Z1
?1
(
y!y ) dx dy ;
+
(x; y )e?i(x!x+y!y ) dx dy:
can not be the zero function, hence, we are guaranteed that there exist no neighborhood at which the Fourier transform of is zero, since is of compact support. Thus the 2D polynomial must be the zero function, i.e. ak = 0 8k, thus proving that the l + 1 functions are independent. 2 It goes without saying that Theorem 3 can be applied to each of the M ? 1 wavelet functions - p(x). Note also, that the proof depends on the orthogonality of the functions j;n(x), and does not hold for wavelets that are only orthogonal in the lter-bank sense, i.e. that are only lossless. The 1D proofs given in [9, 8], use a completely di erent starting point and are given also for non-orthogonal wavelets, but unfortunately have no simple 2D generalization. The speci c examples given in [22, 21] indicate that choosing the sampling matrix so that Dn = I , results indeed in more regular wavelets. 13
4 Designing wavelets with vanishing moments We now turn to the actual design of orthogonal factorable wavelets. A Rank M wavelet matrix is constructed from, paraunitary building blocks of degree one. With the wavelet matrix at hand, the wavelet lters are reconstructed from their polyphase components; (H 0(z1; z2); H 1(z1 ; z2); : : :; H M ?1(z1; z2)) = (z?k ; z?k ; : : :; z?kM ? )HI (zD )H(M ) ; (32) where ki is the shift vector corresponding to the i'th coset. The general form of the Haar wavelet matrix is ! 1 0 H(M ) = H(cM ) (33) 0 OM ?1 where OM ?1 is an orthogonal (M ? 1) (M ? 1) matrix and H(cM ) is the canonical Mband Haar wavelet matrix [17]. The rst column of the canonical Haar wavelet matrix is p1 (1; 1; : : :; 1)T , hence the scaling lter is given by M 0
1
1
0 1 BB 11 CC H (z ; z ) = (z?k ; z?k ; : : :; z?kM ? )HI (zD ) p1 B B .. CCC : MB @.A 0
1
2
0
1
1
(34)
1 Thus a factorable scaling lter is de ned ! by (m + n) unit vectors, that is by (M ? 1)(m + n) M ? 1 parameters which de ne the Haar wavelet matrix free parameters. The other 2 do not e ect the scaling lter, but only the other wavelet lters. Using any of the three alternatives of Theorem ! ! 2, one can see that for L vanishing L +1 moments we have Nc = (M ? 1) ? 1 constraints, hence in general we need 2 wavelets with at least Nfp = Nc free parameters. If one uses the third alternative of Theorem 2, that is equations (17), it is obvious that only the parameters de ning the scaling lter can be taken into account. The other parameters do not appear in equation (34) and therefore do not appear in the constraint equations (17) either. This decoupling of the Haar wavelet parameters from the vanishing moment constraints is the reason for de ning the polyphase (wavelet) matrix the way we did, instead of the more common way which uses the transpose of the one we use. (Another way of reaching this goal is to use a di erent factorization, in which the characteristic matrix is the rst factor on the left, and not on the right.) These considerations suggest that in order to obtain wavelets with ! L vanishing moments we should L + 1 ? 1. Furthermore, in both cases use wavelets with McMillan degrees (m + n) 2 dealt with, the two-band quincunx and four-band separable sampling, it is obvious from 1D results that each polynomial degree should be larger or equal to L ? 1. Otherwise we could nd 1D projections with L vanishing moments and length less than 2L, contrary to 1D well known results [5]. The rest of the section is devoted to two types of orthogonal 2D wavelets, quincunx and four-band separable sampling wavelets. 14
4.1 Quincunx Wavelets
The quincunx sampling is a 2D, two-band, nonseparable sampling scheme. Detailed analysis of 2D multirate systems with quincunx sampling can be found in Kovacevic's doctoral thesis [3]. The sampling and Haar wavelet matrices we use are (there are many other possibilties):
!
1 1 ; 1 ?1
D=
H = p12 (2)
!
1 1 ; 1 ?1
and choose ki = f(0; 0)T ; (1; 0)T g as the coset shift vectors. With this choice of the sampling matrix we have D2 = 2I , and it has been shown to result in wavelets with higher regularity than with other sampling matrices[3, 22, 21]. In the quincunx case the McMillan degrees (m; n) of the wavelet matrix are equal to the polynomial degrees of the elements of the matrix, so for the rest of the subsection we simply use the term degree. Similar to the two-band 1D case, the wavelet lter is fully determined by the scaling lter. For a scaling lter given by: H 0(z1; z2) = H00(z1z2 ; z1z2?1 ) + z1?1 H10(z1z2 ; z1z2?1 ); the wavelet lter is:
H 1(z1 ; z2) = (z1 z2)?m (z1z2?1 )?n H10 (z1?1z2?1 ; z1?1z2 ) ? z1?1(z1z2)?m(z1z2?1)?n H00(z1?1z2?1; z1?1z2): In the two-band factorable case the vectors Vk in equation (9) are 2D unit vectors and thus can be parameterized by one parameter, V (k ) = (cos(k ); sin(k ))T , where k 2 [0; 2 ): The general form for the vanishing moment constraints on the angles k that de ne the wavelet system is developed in Appendix A. We use the third alternative in Theorem 2, i.e. , we require the frequency response of the scaling lter to have a zero ! of order L at the L + 1 ? 1 constraints aliasing point (z1 = ?1; z2 = ?1). We get the following Nc = 2
p+l 0 0 = @ Hp (z1l; z2) @z @z X 1 2 m;nz =z =?1 = ?2?k;q (p; l)Sk;q + p?m;n k;q (p ? 1; l)(Sk;q + Ck;q ) 1
k q even
X
( + )
+
2
k q odd
( + )
m;n 2?m;n k;q (p; l)Ck;q + p?k;q (p ? 1; l)(Sk;q ? Ck;q );
(35)
with (p + l) = 1; 2; : : :; (L ? 1). The coecients ?m;n k;q (p; l) were computed and tabulated, +n . As mentioned earlier, and Sk;q and Ck;q are trigonometric functions of the angles fgmk=1 in order to achieve L vanishing moments one should use in general wavelets of degrees (m + n) Nc . Since the wavelets with vanishing moments are obtained by solving a set of non-linear equations, in general we will have more than one solution with minimal degrees, for a given number of vanishing moments. This is true for the 1D case as well, but contrary 15
to the 1D wavelets with vanishing moments dealt by Daubechies, here di erent solutions may have completely di erent spectral characteristics. We would like to mention the special case of solutions that are actually 1D lters.
Lemma 2 A degree L?1 1D wavelet lter with L vanishing moments is likewise a factorable quincunx wavelet lter with degrees (L ? 1; L ? 1) and L vanishing moments. This can easily be demonstrated by factorizing each degree-1 factor of the 1D polyphase matrix as
I + (z12 ? 1)ViViT = I + (z1z2 ? 1)ViViT I + (z1z2?1 ? 1)ViViT :
These 1D quincunx wavelets were suggested by Cohen and Daubechies[22]. Although the wavelet lters are 1D, due to the non-separable sub-sampling the resulting wavelet functions are 2D non-separable functions. Thus, the above arguments concerning the minimal degrees needed in order to obtain L vanishing moments are true in general, but for the special case of 1D lters we can satisfy the constraints, using wavelets of degrees (L ? 1; L ? 1) even though 2(L ? 1) < Nc for all L > 2. We now give a few examples, but rst note that:
The vectors Vi appear in H(z ; z ) (see equation (9)) only in the combination ViViT , 1
hence, since
2
V (i )V (i )T = V ( + i )V ( + i)T we can limit the search to i 2 [? 2 ; 2 ). for every solution fig there exists another solution f( ? i)g, that will give rise to a scaling lter (function) which is the re ection of the original lter (function) through its center, i.e. , ~h0 (n1 ; n2) = h0 (N1 ? n1 ; N2 ? n2 ).
Given a solution h (n ; n ) with degrees (m; n) there exists another solution with degrees (n; m) in which the variable z is interchanged with z in the polyphase matrix. This will cause the transformation hi (n ; n ) ! hi (n ; ?n ). Hence the following examples are restricted to degrees (m; n) where m n. Example 1.1: The smallest quincunx wavelet lter is the Haar wavelet, its degrees are (m; n) = (0; 0), thus the corresponding wavelet matrix is the Haar wavelet matrix H . 0
1
2
1
2
1
2
1
2
(2)
The lters are exactly the 1D Haar lters, i.e. h0(n1 ; n2) := p1 (1 ; 1) and h1(n1; n2) := p1 (1 ; ?1); 2 2 and only the rst moment vanishes. The absolute value of the frequency response of the scaling lter, and the scaling function in this case are plotted in Fig. 1.a and 1.b respectively. Example 1.2: We next deal with wavelets with L = 2 vanishing moments, in this case the minimal degrees are (1; 1). The two constraint equations are: 1.
0 = 1 ? cos(21 ) ? cos(22) 16
2. 0 = cos(21 ) ? cos(22): We have two di erent solutions (and their re ections), 1. 21 = 22 = 3 : which in accordance with Lemma 2 is the 1D, 4-tap lter of Daubechies with two vanishing moments.
1 0 0 0 0 0 B p p p p C h (n ; n ) := p1 B 1+ 3 3+ 3 3? 3 1? 3 C A @ 4 2 0
1
2
0
0
0
0
2. 21 = ?22 = 3 : which produces the scaling lter
1 0 p p 3 3 + 3 0 0 3 + 3 B p p p p C: h (n ; n ) := p1 B ? 1 ? 3 3 + p3 3 ? p3 ?1 + 3 C A @ 8 2 0 0 3? 3 3?3 3 0
1
2
The frequency response of the scaling lter and the scaling function for the rst solution are shown in Fig 1.c and 1.d, and the corresponding functions for the second solution are depicted in Fig. 1.e and 1.f. For L > 2 the constraint equations can not be solved analytically and a numerical program given in Numerical Recipes [24] was used for the solution of the non-linear equations. ! m + n In addition to the (m + n) free parameters de ning the wavelet matrix, there are
n
alternatives in which the factors of degree-1 can be arranged in equation (9). Each alternative results in a di erent set of constraints, thus introducing another discrete parameter into the problem. It turns out that a solution does not exist for all possible alternatives. For example, analyzing the equations for wavelet lters with degrees (3; 2) we found that one can not obtain three vanishing moments in the two cases in which all the factors of degree-1 in one variable precede the factors in the other variable. Example 1.3: The minimal degrees for obtaining three vanishing moments (not including 1D lters) are (3; 2). The functions plotted in Fig. 2.a and 2.b correspond to a solution obtained by choosing the combination (z1; z1; z2; z1; z2) and the ve parameters as
f2gi = f?0:896794; ?2:934522; 0:782106; 0:546171; 1:781789g: 5 =1
In Figures 2.c and 2.d we show another solution, this solution was obtained using the order (z1 ; z1; z2; z2; z1), and the parameters
f2gi = f0:808218; 1:258868; 0:782106; 1:781789; ?2:091721g: 5 =1
Although both solutions have three vanishing moments there is a considerable di erence between them. The rst solutions is much closer to an ideal quincunx low-pass lter, with 75% of its energy concentrated in the low-pass band, as opposed to the second solution with 60%. The di erence becomes also apparent by looking at the scaling functions. The 17
regularity measure, introduced in the previous section, gave the value = 0:6 0:1 for the rst solution, and = 0:05 0:1 for the second. We calculated the lower bound of for the second example. Using the method suggested by Villemoes[10], we were not able to show that the function is even continuous ( > 0), when applying the method for resolution levels j 14. The 1D solution, corresponding to the Daubechies wavelet with three vanishing moments, obtained the highest value of regularity with = 1:0 0:1, which agrees with the 1D bound for the Holder exponent of = 1:0878. For the sake of completeness we searched the entire space of wavelet matrices with degrees (3; 2) and three vanishing moments, we did so using the complete charicterazation of 2D PUP matrices given in [12]. To obtain compactly supported wavelets with three vanishing moments, one must solve 15 equations for 15 variables. These equations were solved numerically, using the same procedure mentioned above. All the solutions found coincide with solutions that were obtained using factorable wavelets, but it is feasible that other solutions, that were not found numerically, do exist. Finally we give two examples of quincunx scaling lters and functions with four vanishing moments. The Frequency response of the scaling lter and the scaling function are plotted in Fig. 3.a and 3.b for one example and in Fig. 3.c and 3.d for the other. The rst example has a good frequency response with 80% of the energy in the low-pass band, but is very spiky, = 0:5 0:1. While the second example is more regular, with = 0:9 0:1, but with a slightly inferior frequency response.
4.2 Four-band, separable sampling wavelets
This sampling scheme is the most commonly used in all applications of the 2D DWT , since it is the simplest ! generalization of the two-band 1D sampling. The sampling matrix is 2 0 D= and the four-band Haar wavelet matrix has the general form; 0 2
0 1 p 1 ? 3 0 0 q B CC p ! B 3 ? 2 0 1 1 = B C 1 0 1 q B C p ; p C H = 2B 2 0 O 1 1 = 3 ? B C @ p q p A 2 3
(4)
1 1= 3
2 3 2 3
3
(36)
2
O3 being a rank-3 orthogonal matrix. Thus, we have a three parameter space of Haar wavelet matrices. We chose the four coset shift vectors to be: ki = f(0; 0)T ; (1; 0)T ; (0; 1)T ; (1; 1)T g, hence, the three points of repeated spectra are: ! i = f(; 0)T ; (0; )T ; (; )T g: With this choice the expression for a factorable scaling lter of McMillan degree (m; n) is according to (34):
0 1 1C B m n B C Y H (z ; z ) = 21 (1; z ; z ; z z ) [I + (zk? ? 1)Vk VkT ] B BB 11 CCC : k @ A 0
+
1
2
1
2
1 2
=1
18
2
1
(37)
Once again, in order to obtain L vanishing moments we impose Nc = 3
!
L+1 ?1
!
2 constraints on the 3(m + n) free parameters de ning the scaling lter. Just as we did in the quincunx case, we seek solutions with the minimal ! possible degrees which in general are L + 1 ? 1. Once we have designed the scaling given by Nfp = Nc , that is (m + n) = 2 function by xing the 3(m + n) free parameters we still have three free parameters de ning the Haar wavelet matrix. Note that using 1D wavelets with L vanishing moments to design separable 2D wavelets will result in wavelets with McMillan degrees ! (2L ? 2; 2L ? 2) (see section 2.2). Since in L + 1 ? 1 for all L 6= 6, in the following examples we this case (m + n) = 4(L ? 1) 6= 2 will not obtain separable lters as private solutions. Example 2.1: The Harr wavelet matrix gives rise to the smallest wavelet lters. Each lter has four non-zero coecients, and the corresponding function is compactly supported on [0; 1]2. In this case only the rst moment vanishes and the resulting functions are discontinuous and piecewise constant. Example 2.2: In order to obtain wavelets with L = 2 vanishing moments we must use wavelets with McMillan degrees (1,1), thus the wavelet matrix is given by,
H(z ; z ) = [I + (z? ? 1)V V T ][I + (z? ? 1)V V T ]H : 1
2
1
1
1 1
2
1
2 2
(4)
We have six bilinear constraints, (two derivatives at three aliasing frequencies) on eight unknowns, (two vectors of length 4) de ned by six free parameters. The two unit vectors solving this set of equations are, p p p p V1 = p1 ( 3; 1; 3; 1)T ; V2 = p1 ( 3; 3; 1; 1)T : 2 2 2 2 As mentioned previously, the three parameters de ning the Haar wavelet matrix do not in uence the number of vanishing moments but rather control, in some sense, the orientation of the three wavelet lters/functions. In the following example we used the separable Haar wavelet matrix 1 0 1 1 1 1 C BB 1 ?1 1 ?1 C 1 C B (4) H = 2B B@ 1 1 ?1 ?1 CCA ; 1 ?1 ?1 1 this is one way to guarantee that each wavelet lter has a maxima at one aliasing point and a zero at all others. Once again, for L > 2 only numerical solutions can be obtained. In Fig. 4 an example of wavelet lters and functions with four vanishing moments are depicted. Their McMillan degrees are (5; 4), the lters are of length (12; 10), and the regularity of the scaling function was measured to be = 1:25 0:1. 19
5 Discussion and Conclusions Several methods have been suggested for constructing 2D, orthogonal wavelet lters and functions [1, 22]. These methods are based on some sort of transformation from 1D wavelets. Hence, the resulting wavelets in each case, constitute only a small restricted subset of all possible orthogonal wavelets. Here we use a parameterization of all factorable orthogonal wavelets, which is a much more general subset. The subset of factorable wavelets was shown in [12] to include almost all possible orthogonal wavelets, and contrary to the complete set of orthogonal wavelet, has a simple parameterization. We have shown that this subset includes the previously suggested 2d orthogonal wavelets, such as, separable wavelets and the special (actually 1D) quincunx wavelets [22]. In order to obtain wavelets with several vanishing moments, one has, in general, to solve numerically a set of non-linear equations. Even when focusing our attention only on wavelets with maximal number of vanishing moments for given McMillan degrees, considerable freedom is still left, since the set of non-linear equations have many solutions. Furthermore, usually, satisfying the conditions of having a large number of vanishing moments is not enough, as is illustrated in Fig. 2 and Fig. 3, Extra care should be taken in choosing the best solution for a given 2D application. Theorem 4 and our examples show that, a sucient number of vanishing moments is a necessary, but by no means sucient, condition for the regularity of the wavelet function. The regularity of the original Daubechies orthogonal wavelets (minimal phase wavelets) was found to grow at least linearly with the number of vanishing moments, yet in 2D we were not able to select any special subset of solutions and to give a lower bound for their regularity. If a high degree of regularity is required for given McMillan degrees, one should probably sacri ce some of the vanishing moments and use the extra degrees of freedom to optimize the regularity. For a given order of moments, one can relax only part of the partial derivative constraints (17), and only at some of the aliasing points. Our treatment of the regularity issue is quite preliminary and mainly serves as a motivation for investigating 2D wavelets with vanishing moments. Further work is required in order to establish the exact relation between the number of vanishing moments and regularity, and its dependency on the choice of sampling matrix and coset vectors. As in many other areas of mathematics, most of the diculties encountered when trying to generalize 1D theory to higher dimensions arise already at the transition from 1D to 2D. Once these diculties have been overcome, the design of N-dimensional wavelets will be a simple generalization of the 2D method.
Appendix A In this appendix we will develop the constraint equations for vanishing moments in the case of factorable quincunx wavelets. Each factor of degree-1 is de ned by one parameter i
Wi = (I + (zi? ? 1)ViViT ); 1
20
were zi = z1 or z2 and Vi is a unit vector of length 2. If we use Vi = (cos(i ); sin(i ))T ; we get !! cos(2 ) sin(2 ) 1 i i (A:1) Wi(i) = 2 (zi?1 + 1)I + (zi?1 ? 1) sin(2 ) ? cos(2 ) : i i Next we de ne two matrices, the 2D rotation matrix- R() and the rotation-plus-re ection matrix T () , ! ! cos(2) ? sin(2) cos(2) sin(2) ; T () = ; R() = sin(2) cos(2) sin(2) ? cos(2) and recall the relations: T ( )T ( ) = R( ? ); R( )R( ) = R( + ); R( )T ( ) = T ( + ); T ( )R( ) = T ( ? ): (A:2) From this we see that the 2D polyphase matrix of McMillan degrees (m; n) can be written as : mY +n m X n X m;n(z ; z )T (fg); H(z1; z2) = W(i) = (A:3)
k;q 1 2 k;q k=0 q=0
i=1
where we de ne 4 m;n(z ; z ) =
k;q 1 2
1 m
n
+
(z1?1 + 1)m?k (z1?1 ? 1)k (z2?1 + 1)n?q (z2?1 ? 1)q ; 2 and Tk;q (fg) is the sum of all possible combinations of an ordered product of k T -matrices of z1 factors and q T -matrices of z2 factors. With use of the relations (A.2) one can see that for a given k and q , Tk;q (fg) will include terms of the following type
8 > > > 4< Tk;q (fg) = > > :
! ?Sk;q =4 P R( ? + ? ? ); if (k+q) is even j j j jk q Ck;q ! Sk;q =4 P T ( ? + ? + ); if (k+q) is odd j j j jk q ?Ck;q For example with degrees (m; n) = (2; 2) ,where the rst two factors are in z and the last in z , we have 16H(z ; z ) = (z ? + 1) (z ? + 1) + (z ? ? 1)(z ? + 1) [T ( ) + T ( )] + (z ? + 1) (z ? ? 1)[T ( ) + T ( )] + (z ? ? 1)(z ? ? 1)[R( ? ) + R( ? ) + R( ? ) + R( ? )] + (z ? ? 1) (z ? + 1) [R( ? )] + (z ? + 1) (z ? ? 1) [R( ? )] + (z ? ? 1) (z ? ? 1)[T ( ? + ) + T ( ? + )] + (z ? ? 1)(z ? ? 1) [T ( ? + ) + T ( ? + )] + (z ? ? 1) (z ? ? 1) [R( ? + ? )]: Ck;q Sk;q Ck;q Sk;q
1
2
3
+
1
2
3
+
1
2
1
2
1
1
1
1
1
1
1
1
1
1
2
2
1
2
2
2
2
1
2
1
2
1
2
2
1
2 2 2
2
2
2
2 1
2 1 2
1
2
2
2
1
2
3
4
1
1
2
1
2
2
2
1
2
3
1
1
2
3
4
4
2
3
1
2
3
1
2
4
1
3
4
2
3
4
1
2
21
3
4
2
4
With an explicit expression for H(z1; z2) we can obtain the vanishing moments constraints as constraints on the parameters . First we construct the scaling lter from the polyphase matrix ! 1 1 ; (A:4) H 0(z1 ; z2) = (1; z1)H(z1z2 ; z1z2?1 ) p 2 1 and then we write the partial derivatives of the scaling function
p @ p lH (z ; z ) 2 @z p @z l z +
0
1
m X n h X k=0 q=0
1
2
2
z ?1
=
1= 2=
i
m;n ?m;n k;q (p; l)(Ck;q(1 1) ? Sk;q (1 1)) + p?k;q (p ? 1; l)(Sk;q Ck;q ) ;
where the upper sign should be used when (k + q ) is even and the lower sign otherwise, and we used the notation m;n ?1 4 @ p+l k;q (z1 z2 ; z1z2 ) m;n : ?k;q (p; l) = p l
@z1 @z2
z1 =z2 =?1
Finally we arrive at the constraints for L vanishing moment. ! By equating all partial derivaL +1 tives (p + l) = 1; 2; : : :; (L ? 1) to zero we get ? 1 constraints 2
i X h m;n ?2?k;q (p; l)Sk;q + p?m;n ( p ? 1 ; l )( S + C ) k;q k;q k;q k q even i X h m;n + 2?k;q (p; l)Ck;q + p?m;n k;q (p ? 1; l)(Sk;q ? Ck;q ) :
0 =
( + )
k q odd
(A.5)
( + )
Acknowledgments
We would like to thank N. Cohen, M. Zibulski of the Department of Electrical Engineering, and Y. Binyamini of the Department of Mathematics, Technion, for their help. We also wish to thank L. Villemoes for his help in applying his methods to our examples.
References [1] S. Mallat, \A theory of multiresolution signal decomposition: The wavelet representation," IEEE Trans. PAMI, vol. 11, pp. 674{693, 1989. [2] T. Chen and P. P. Vaidyanthan, \Mulidimensional multirate lters and lter banks derived from one-dimensional lters," IEEE trans. on SP, vol. 41, May 1993. [3] J. Kovacevic, Filter banks and Wavelets: Extentions and Applications. PhD thesis, Department of Electrical Engineering, Columbia University, 1991. 22
[4] C. Lien and C. Huang, \The design of 2-d nonseparable directional oerfect reconstruction lter banks," Multidimensional Systems and Signal Processing, vol. 5, p. 289, 1994. [5] I. Daubechies, \Orthonormal bases of compactly supported wavelets," Comm. Pure and Appl. Math, vol. 41, pp. 909{996, 1988. [6] O. Rioul, \On the choice of wavelet lters for still image compression," in ICASSP, vol. V, p. 550, 1993. [7] I. Daubechies, Ten lectures on wavelets. CBMS/NSF series in applied math., SIAM, 1992. [8] I. Daubechies and J. C. Lagarias, \Two-scale di erence equations II. Local regularity, in nite products of matrices and fractals," SIAM J. Math. Anal., vol. 23, no. 4, pp. 1031{1079, 1992. [9] O. Rioul, \Simple regularity criteria for subdivision schemes," SIAM J. Math. Anal., vol. 23, no. 6, pp. 1544{1576, 1992. [10] L. F. Villemoes, \Continuity of nonseparable quincunx wavelets," Appl. Comput. Harmon. Anal., vol. 1, pp. 180{187, 1994. [11] G. Karlsson and M. vetterli, \Theory of two dimensional multirate lter banks," IEEE trans. on ASSP, vol. 38, June 1990. [12] S. Venkataraman and B. C. Levy, \State space representations of two dimensional FIR lossless transfer matrices," IEEE Trans. Circuits and Syst., vol. 41, p. 117, Feb. 1994. [13] S. Venkataraman and B. C. Levy, \A comparison of design methods for 2-D FIR orthogonal perfect reconstruction lter banks." submitted to IEEE Trans. Circuits and Syst., 1994. [14] J. Kovacevic and M. Vetterli, \Nonseparable two and three-dimensional wavelets," IEEE Trans. Signal Process., vol. 43, May 1995. [15] E. Viscito and J. P. Allebach, \The analysis and design of multidimensional FIR perfect reconstruction lter banks for arbitrary sampling lattices," IEEE Trans. Circuits and Syst., vol. 38, pp. 29{41, Jan. 1991. [16] W. M. Lawton and H. L. Resniko , \Multidimensional wavelet bases," tech. rep., Aware, Aware Inc. Cambridge MA, 1991. [17] P. N. Heller, H. L. Resniko , and R. O. Wells, \Wavelet matrices and the representation of discrete functions," in Wavelets: a tutorial in theory and applications (C. K. Chui, ed.), Academic press, 1992. [18] P. P. Vaidyanthan, T. Q. Nguyen, Z. Doganata, and T. Saramaki, \Improved technique for design of perfect reconstruction FIR QMF banks with lossless polyphas matrices," IEEE Trans. on ASSP, vol. 37, July 1989. 23
[19] D. Stanhill and Y. Y. Zeevi, \Two-dimensional orthogonal wavelets, with linear phase," Tech. Rep. CC PUB#111, Dep. of Electrical Eng. , Technion, Haifa, ISRAEL, 1995. [20] P. Ste en et al., \Theory of regular m-band wavelet bases," Tech. Rep. TR-93-02, Computational Math. Lab., Rice University, Houston, TX, 1993. [21] K. Grochenig and W. R. Madych, \Multiresolution analysis, Haar bases, and selfsimilar tilings of Rn ," IEEE Trans. Inform Theory, vol. 38, March 1992. [22] A. Cohen and I. Daubechies, \Non-separable bidimensional wavelet bases," Rev. Mat. Iberoamericana, vol. 9, no. 1, pp. 51{137, 1993. [23] G. G. Walter, Wavelets and other orthogonal systems with applications. CRC Press, 1994. [24] W. H. Press et al., Numerical Recipes in C. Cambridge University Press, 1992.
FIGURE CAPTIONS Fig. 1: Quincunx scaling lters and functions. The frequency response of the scaling
lter, and the scaling function of: (a) and (b) the Haar quincunx wavelet (Example 1.1), (c) and (d) wavelet with two vanishing moments, 1D solution (Example 1.2), (e) and (f) wavelet with two vanishing moments, second solution (Example 1.2).
Fig. 2: Quincunx scaling lters and functions with three vanishing moments. The
frequency response of the scaling lter, and the scaling function of the rst ((a) and (b)) and second ((c) and (d)) solutions of Example 1.3. Both scaling functions show the eighth order approximation, using the sub-division scheme.
Fig. 3: Quincunx scaling lters and functions with four vanishing moments. The frequency response of the scaling lter, and the scaling function of the two examples mentioned in the text.
Fig. 4: Four-band, separable sampling wavelet with four vanishing moments. Figures (a)-(d) show the magnitude response of the four wavelet lters, and (e)-(h) show the wavelet functions.
24
Abstract We investigate a very general subset of two-dimensional, orthogonal, compactly supported wavelets. This subset includes all the wavelets with a corresponding wavelet (polyphase) matrix that can be factored as a product of factors of degree-1 in one variable. In this paper we consider in particular wavelets with vanishing moments. The number of vanishing moments that can be achieved increases with the increase of the McMillan degrees of the wavelet matrix. We design wavelets with the maximal number of vanishing moments for given McMillan degrees, by solving a set of nonlinear constraints on the free parameters de ning the wavelet matrix, and discuss their relation to regular, smooth wavelets. Design examples are given for two fundamental sampling schemes; the quincunx and the four-band separable sampling. Their relation to the, well known, one-dimensional Daubechies wavelets with vanishing moments is discussed.
This research was supported in part by the Ollendorf Center of the Department of Electrical Engineering, by the Israel Science Foundation and by the Fund for Promotion of Research at the Technion.
1
1 Introduction The M-band 2D discrete wavelet transform (DWT) is de ned by M functions, one scaling function - 0(x) and (M-1) wavelet functions - p(x), and a sampling matrix - D. For wavelets associated with multiresolution analysis the DWT can equally be speci ed by M wavelet lters - H p(!), instead of the M wavelet functions. In this work we restrict the discussion to compactly supported functions that de ne an orthogonal DWT. The orthogonal compactly supported wavelet lters are a special case of lossless FIR multirate lter banks. Since considerable work has been carried out on lossless lter banks, we adopt, in part, the notation and the terminology used in this eld. Speci cally, we limit most of the analysis to the polyphase matrix de ned by the wavelet lters (wavelet matrix) and not to the wavelet functions. Most current applications of the 2D DWT use 1D wavelet functions and lters to construct 2D wavelets, usually separable wavelets[1]. Using separable functions preserves orthogonality but imposes a severe limitation on the 2D wavelets. For example the commonly used four-band lters can not have linear phase, except for the Haar wavelet, and their orientation can not be adjusted for given speci cations. Various other methods were proposed for the design of 2D perfect reconstruction lter banks [2, 3, 4], which do not require orthogonality (paraunitarity) or nite support. The restrictions imposed on the functions/ lters still leave considerable freedom in choosing the wavelet functions/ lters. Motivated by the work of Daubchies in one dimension[5], we examine the design of 2D, orthogonal, compactly supported wavelets with maximal number of vanishing moments. In the 1D case these have become the standard wavelets, and they can be considered suboptimal in a large class of applications, such as the representation, processing and compression of 1D signals. For detailed experimental results in image compression using separable wavelets, showing the optimality of wavelets with vanishing moments, see [6]. The generalization of the 1D notion of vanishing moments to higher dimensions is straightforward, but its nature and importance is still to be determined. In 1D 2-band wavelets, the number of vanishing moments is closely related to the atness of the frequency response of the wavelet lters at the aliasing points. When dealing with more than two bands, even in 1D, the M functions are no longer uniquely de ned by each other, and the vanishing moment constraints ensure atness of the scaling lter but not of the other wavelet lters. The regularity of the wavelet functions, is also closely related to the number of vanishing moments in the 1D case[7, 8, 9]. Yet, in more than 1D no such relation has been rigorously shown, although it was suggested to hold by Daubechies and Lagarias[8]. The work done by Villemoes on the quincunx sampling scheme[10], introduces a di erent approach to the general 2D regularity measurement, and is used here to obtain lower bounds of the Holder exponent. The method used by Daubechies[5] in 1D can not be generalized to construct orthogonal wavelets with N vanishing moments in more than one dimension. The factorization of the frequency response of the scaling lter as 1 + exp(?i!) N 0 L(!); H (!) = 2 has no strait-forward generalization to more than 1D, since the fundamental theorem of 2
algebra holds only for 1D. We therefore construct 2D wavelets by rst parameterizing the space of orthogonal wavelets with compact support, and then impose a set of constraints on the parameters in order to achieve the desired number of vanishing moments. The shortcome of this method is that no complete parameterization of the space of 2D, orthogonal, compactly supported wavelets exists. In most of this work we limit the design to wavelets with wavelet matrices that are factorable as a product of 1D degree-1 factors. Such wavelets were proposed by Karlsson and Vetterli[11]. The issue of parameterization of 2D lossless FIR transfer matrices has been investigated by Venkataraman and Levy[12, 13], they have shown that almost all lossless FIR transfer matrices (wavelet matrices in our case) are factorable. Towards the completion of this manuscript our attention was drawn to a recent work by Kovacevic and Vetterli[14]. They use a construction method similar to the one proposed here, and give a few design examples. The rest of the paper is organized as follows. The rst section is devoted to the description of 2D orthogonal wavelets and lter banks in terms of their polyphase matrix. Two special subclasses are discussed, factorable and separable wavelets and lter banks. In the next section we give the conditions for obtaining wavelets with L vanishing moments, and discuss their relation to regular, smooth wavelets. Speci c design results and examples are given in the last section; rst for quincunx wavelets and then for four-band, separable sampling wavelets. A few comments on notation. The term degree refers to the polynomial degree if not otherwise speci ed. We use boldface for two-dimensional entities like x; n. We adopt the notation used in the multidimensional multirate literature; for two vectors n = (n1 ; n2)T 4 n n and z = (z1 ; z2)T we have zn = z1 z2 , and for a 2 2 matrix D we de ne, zD =4 (z1D ; z2D ; ; z1D ; z2D ; ). Capital letters are used for functions given in the z and frequency domain, such as H (z) and H (!). Finally, boldface is used also for matrix valued functions. 1
11
21
12
2
22
2 2D Wavelets and multirate lter-banks When dealing with multi-dimensional wavelets and multirate systems, the change in resolution and sampling rate is given by a matrix -D; a 2 2 matrix in the 2D case as compared with a single factor in the 1D case. For a detailed analysis of 2D sub-sampling and multirate systems see Viscito and Allebach[15], and for a more mathematical background in multidimensional wavelets see Lawton and Resniko [16]. Here we quote a few facts and de nitions relevant to the rest of this paper. For a discrete signal given on a lattice , the sub-sampled signal is de ned on the lattice D = fDn : n 2 g. The density of grid points on D is smaller by a factor of M = j det(D)j than the density on . A coset of a sub-lattice is the set of grid points on the lattice which are also on a shifted version of the sub-lattice, a coset is referred to by a vector ki , that de nes the 2D shift. There are M distinct cosets which cover the whole lattice. The points !i = 2D?T ki in the frequency plane are called points of repeated spectra or \aliasing points". We now de ne 2D wavelets. The discrete 2D wavelet transform (DWT) of a function f 2 L2 (IR2) is de ned as
Z
4 f (x) p (x) dx; (W f )pj;n = j;n
3
(1)
where
4p j p j p 2 j;n(x) = M (D x ? n) : p = 1; : : :M ? 1; j 2 ZZ ; n 2 ZZ : (M ? 1) wavelet functions - p(x) are related to the scaling function - 0(x) through
The the two-scale re nement equations:
p (x) = pM X hp (n)
(Dx ? n); (2) n where hp (n) is a 2D sequence. In this work we focus our e ort on real, compactly supported wavelets, hence the sequences hp (n) are real and of nite length. The M sequences hp (n) and their Fourier transform H p(!) can be viewed as the impulse and frequency responses of 2D FIR lters, and are called wavelet lters. H 0(!) the scaling lter is a low-pass lter and the rest are high-pass or band-pass lters. It is by means of these wavelet lters, that we obtain, the now well established, connection between wavelet theory and multirate lter banks. To assist the analysis of multirate lter banks we de ne two matrices, the polyphase and modulation matrices. In multirate systems each lter H p(!) is decomposed into M polyphase components, one for each coset
H p(!) =
MX ?1 j =0
0
e?i!kj Hjp(DT !);
(3)
where each polyphase component is given by,
Hjp(!) =4
X
hp (Dn + kj ) e?i!n : (4) n The polyphase matrix of a maximally decimated lter bank is a M M matrix, X (5) H(!) = h(n) e?i!n; n with one column for each lter and one row for each polyphase component, that is the p'th column is T p H0 (!); H1p(!); HMp ?1(!) (Note that the polyphase matrix is usually de ned as the transpose of the one de ned above, we use our de nition for future simpli cation.) The modulation matrix - HM (!) is closely related to the polyphase matrix. In this case the p'th column is built from M shifted versions of the lter H p(!), one for each aliasing frequency, that is the p`th column is T p p M H (!); H p(! + 2D?T k1 ); H p (! + 2D?T kM ?1) :
Claim 1 A maximally decimated lter bank is lossless if and only if it obeys one of the
following equivalent conditions : 1.
P hp(n)hq(Dm + n) = : p;q m; n 0
4
2. Hy(!)H(!) = I or
P hy(n + k)h(n) = I : k; n 0
3. HyM (!)HM (!) = I: Claim 2 A bank of M , 2D, FIR lters de nes a compactly supported, orthogonal 2D wavelet system if and only if : The lter bank is lossless, p Pn hp(n) = H p(! = 0) = p;0 M: The term orthogonal applies only to the wavelet lters and not to the wavelet functions, which in general are not orthogonal but constitute a tight frame with frame bounds A = B = 1. In the z-domain the polyphase matrix of FIR (causal) lters is a 2D polynomial matrix in the variables z1?1 and z2?1 . The unitary condition of the polyphase matrix is now HT (z?1)H(z) = H(z)HT (z?1) = I: (6) Polyphase matrices that obey the above condition are called paraunitary, i.e. unitary at all points jzi j = 1. We use the term wavelet matrix of rank M for the polyphase matrix of an orthogonal M -band wavelet system. Since the characteristic matrix H0 of a paraunitary polynomial matrix (PUP matrix), de ned as H0 = H(z1 = 1; z2 = 1), is a constant orthogonal matrix, any PUP matrix can always be decomposed as, H(z1; z2) = HI (z1; z2)H0; with HI (z1 = 1; z2 = 1) = I: For wavelet matrices the characteristic matrix must have a special form so as to satisfy the constraints in claim 2. These matrices are called Haar wavelet matrices[17], and are a generalization of the, well known, two-band Haar matrix. The notation - H(M ) is used for a Haar wavelet matrix of rank M . Throughout the rest of the paper we will work in the z - domain, except in section 4, where the results are better understood in the frequency domain. Most results will be stated for general PUP matrices but the design examples will be limited to wavelet matrices. The reader unfamiliar with 2D polynomial matrices, and notions such as the McMillan degrees should refer to [12], for details. We quote one result from [12] for reference and further consideration. Claim 3 A polynomial matrix is paraunitary and of McMillan degrees (m,n) if and only if; (7) det H(z1 ; z2) = z1?m z2?n :
2.1 Factorable wavelets
Denition 1 A PUP matrix of McMillan degrees (m,n) in z and z is called factorable 1
2
if it can be factored into a product of McMillan degree-1, 1D matrices, m of them in z1 and n in z2 , that is: mY +n H(z1; z2) = Hk (zk ); (8) k=1
where zk = z1 or z2 , and Hk are 1D PUP matrices of McMillan degree-1.
5
Using 1D results[18] we can write any factorable PUP M M matrix of McMillan degrees (m; n) as: mY +n (9) H(z1; z2) = [I + (zk?1 ? 1)VkVkT ]H0; k=1
where H0 is the characteristic matrix of H, Vk are unit vectors of size M , and as above, zk = z1 or z2. As mentioned earlier, for wavelet matrices H0 must satisfy an extra constraint of regularity, and is denoted by H(M ) . The number of free parameters in the above factorization (9) is ! ! M M ? 1 Nfp = + (M ? 1)(m + n); or, for wavelets + (M ? 1)(m + n); 2 2 where the binomial part comes from the characteristic (Haar wavelet) matrix and the rest from the (m + n) unit vectors of length M . In [12] it is shown that this is equal to the number of degrees of freedom of all PUP matrices of rank M and McMillan degrees (m,n), hence the set of factorable PUP matrices is not a small subset of all PUP matrices. Furthermore, in the case that one of the polynomial degrees is equal to one, all the PUP matrices are factorable[2]. In our complementary study [19] we show that a PUP matrix is factorable in the McMillan degree sense, as de ned above, if and only if it is factorable in the polynomial sense. The polynomial factorization is used for the factorization of linear phase PUP matrices.
2.2 Separable Wavelets
In most applications of the 2D DWT, separable wavelets are used. In this section it is shown how separable wavelets t into the general framework of 2D wavelets. Denition 2 We call a MN -band, 2D lter bank separable if it can be separated into a direct (Kronecker) tensor product of 1D M-band and N-band lter banks, that is (H 0(z1; z2); H 1(z1 ; z2); : : :; H MN ?1(z1 ; z2)) = (H~ 0(z1)H 0(z2); : : :; H~ 0(z1 )H N ?1(z2); H~ 1(z1 )H 0(z2); : : :; H~ M ?1(z1 )H N ?1(z2 )) = (H~ 0(z1); H~ 1(z1 ); : : :; H~ M ?1(z1 )) (H 0(z2); H 1(z2 ); : : :; H N ?1(z2)): ?1 p(z2 )gpN=0?1 are M-band and N-band 1D lter banks respectively. where fH~ p(z1)gM p=0 and fH Separable lter banks are a special case of 2D lter banks with separable sampling, i.e. with a diagonal sampling matrix. A straight forward way of de ning the polyphase matrix for separable lter banks is by choosing the coset shift vectors as: ( ! ! ! ! ! !) 0 0 0 1 1 M ? 1 ki = ; ; : : :; ; ; : : :; ; : : :; : 0 1 N ?1 0 N ?1 N?1 Hence the MN lters can be reconstructed from their polyphase components, (H 0(z); H 1(z); : : :H MN ?1(z)) = (zk ; zk ; : : : zkMN ? )H(zD ) = (1; z2; : : :z2N ?1 ; z1; : : :z1M ?1 z2N ?1 )H(z1M z2N ): 0
6
1
1
With this choice of coset shift vectors we have the following lemma. Lemma 1 A 2D MN -band PUP lter bank system is separable if and only if its polyphase matrix (with the above de ned coset shift vectors) is given by: H(z1; z2) = H~ (z1) H (z2); ~ (z1 ) and H (z2) are 1D M-band and N-band PUP matrices respectively. where H Proof: using the identity [A B ][C D] = [AC ] [BD]: (10) we have (1; z2; : : :z2N ?1 ; z1; : : :z1M ?1 z2N ?1 )H(z1M ; z2N ) = (H 0(z1; z2); : : :; H MN ?1(z1; z2)) = (H~ 0(z1); : : : H~ M (z1)) (H 0(z2 ); : : : H N (z2 )) ~ (z1M )] [(1; z2; : : :z2N ?1 )H (z2N )] = [(1; z1; : : :z1M ?1)H (z2N )] ~ (z1M ) H = [(1; z1; : : :z1M ?1) (1; z2; : : :z2N ?1 )][H ~ (z1M ) H (z2N )]; (11) = (1; z2; : : :z2N ?1 ; z1; : : :z1M ?1 z2N ?1)[H ~ (z1 )
which proves that H(z1; z2) is separable only if it can be decomposed as H(z1 ; z2) = H H(z2). The converse can be obtained by reversing the order of the equalities in (11). Using (10) once more, we get H(z1; z2) = H~ (z1) H (z2) ~ I (z1) H I (z2 )][H(0M ) H0(N )]: (12) = [H Hence, separable orthogonal MN -band lter banks have a separable characteristic matrix. We still have to show that paraunitarity is preserved. If both 1D matrices are paraunitary we have, H(z1; z2)H(z1?1; z2?1)T = H~ (z1) H (z2) H~ (z1?1) H (z2?1) T ~ ~ T ?1 T ?1 (z1 )H (z1 ) H(z2 )H (z2 ) = H = I (M ) I (N ) = I (M +N ): The converse is also true (once we x the normalization) since, A B = I (M +N ) ) A = diag(a; a; : : :; a) and B = diag(1=a; 1=a; : : :; 1=a); for some a 6= 0.
2
It is important to note that due to the fact that ~ N ? M ~ H ) = [det(H ~ I (N ))][det(I (M ) (H )] = det(H det(H ) det(H) ; from claim 3 we see that the direct tensor product of two 1D, PUP matrices of McMillan degree m and n, respectively, results in a 2D PUP matrix of McMillan degrees (Nm; Mn). 7
Theorem 1 Every 2D, MN -band, separable PUP matrix constructed from a M -band 1D PUP matrix with McMillan degree m, and a N -band 1D PUP matrix with McMillan degree n, is factorable and of McMillan degrees (Nm; Mn) Proof: Using identity (10) we have; 3 "Y # 2Y m n HI (z ; z ) = [I M + (z ? ? 1)ViViT ] 4 [I N + (z ? ? 1)Uj UjT ]5 i j 3 "Y #2 m m Y = [I M + (z ? ? 1)ViViT ] I N 4I M [I N + (z ? ? 1)Uj UjT ]5 i j n m i i Y h MN Y h MN I + (z ? ? 1)(I M Uj UjT ) = I + (z ? ? 1)(ViViT I N ) j i n m i i Y h MN Y h MN I + (z ? ? 1)Qj ; I + (z ? ? 1)Pi = 1
(
2
)
1
1
(
=1
)
2
1
=1
(
)
1
(
1
)
(
)
(
=1
)
2
1
=1
(
)
1
1
(
(
)
)
2
1
(
)
=1
=1
(
)
1
i=1
1
(
)
2
j =1
1
where Vi and Uj are unit vectors of length M and N respectively, and Pi and Qj are matrices, the nature of which is yet to be determined. De ning the vector k to be a vector of zeros except for the k'th entry which is one, we can write
Pi = ViViT I (N ) = Vi ViT
=
N X
N X k=1
k Tk
(Vi k )(Vi k )T :
k=1
(13)
Since the vectors (Vi k ) ; k = 1; : : :N are unit vectors of length MN and are orthogonal to each other, Pi = ViViT I (N ) is a rank-N symmetric projection matrix from IRMN . In a similar manner we have
Qj = I
M) U
(
T j Uj
=
M X
k=1
(k Uj )(k Uj )T ;
(14)
which again can be seen to be a rank-M symmetric projection matrix from IRMN . Summing up, we have Nm factors of McMillan degree-1 in z1 and Mn factors in z2 , thus completing the proof. It goes with out saying that the polynomial degrees of a separable PUP matrix are equal to the polynomial degree of its two 1D factors in each variable respectively. 2
3 Vanishing Moments and Regularity 1D wavelets with vanishing moments have been used extensively in the context of 1D and 2D applications of the DWT . The requirement of maximal number of vanishing moments was found to be appropriate for a variety of applications, and gave close to optimal results in 8
most cases. This success in 1D motivates us to study 2D wavelets with vanishing moments, and to apply them in image coding and representation, as well as in other 2D applications. As in the 1D case[20] we can characterize the quality of vanishing moments (or Lregularity) in several analogous manners. Theorem 2 . The following three conditions are equivalent, and can serve as a de nition of an orthogonal M-band wavelet system with L vanishing moments. 1. All moments up to order (L - 1) of the wavelet lters vanish, that is
k+q H p(!x ; !y ) X pk;q =4 hp (n)nq = (?i)k+q @ @! = 0; q xk @!y n !=0
(15)
for all p = 1; : : : (M ? 1) and k + q = 0; 1; : : :L ? 1, where nq denotes (m; n)(k;q) = mk nq . 2. All moments up to order (L - 1) of the wavelet functions vanish, that is: Z k+q ^p(!x ; !y ) @ 4 p mk;q = xq p(x) dx = (?i)k+q = 0; k q
@!x @!y
(16)
!=0
for p; k; q as above and again xq denotes xk y q .
3. The frequency response of the scaling lter has a zero of order L at the (M ? 1) aliasing frequencies, i.e. @ k+q H 0(!x ; !y ) (17) k q = 0;
@!x@!y
!=!i
for all (M ? 1) aliasing frequencies ! i and (k + q ) = 0; 1; : : :L ? 1.
Proof: 1 ) 2 : using equation (2) we can derive the recursion equation p X p Z q p mk;q = M h (n) x (Dx ? n) dx n p X p Z ? = M h (n) (D ( + n))q () dM : 0
1
n Writing (D?1 ( + n))q explicitly and using the binomial expansion we arrive at
mpk;q
! !X q q?i?j k X X k q i+j k+X 1 i q ? i = p ab i M i=0 j=0 ps;t m0(i+j?s);(k+q?i?j?t) ;
j
s=0
(18)
0
t=0
i+j s
!
!
k+q?i?j t
(19)
where D?1 = (a; b), and i = (i; j )T . Hence if pk;q = 0; 8(k + q ) < L on the right hand side of (19) we get mpk;q = 0 on the left hand side. 2 ) 1 : We use the recursion relation (19) and prove by way of induction. a) k = q = 0: p0;0 = 0 must be true for all wavelets (claim 2). 9
b) k + q = l: Given mpk;q = 0 8 (k + q ) l, we assume that pk;q = 0 8 (k + q ) < l and will prove that pk;q = 0 8 (k + q ) = l. Using the above assumption, in equation (19) we are left only with the terms in which s + t = l, that is, for k + q = l
! q k X X k = aibq?i i i j 0 Xl @X i q?i = ab
0 = mpk;q
=0 =0
i j
( + )=0
i?j
!
q p m0 j i+j;k+q?i?j 0;0 ! !1 k q A p )m00;0; i +j;l?(i+j ) i j
since m00;0 6= 0, for each p, we are left with l + 1, linear homogeneous equations for l + 1 unknowns - pi+j;l?(i+j ) . Since D is not singular the linear system is non-singular and we arrive at pi;l?i = 0 8i l. 3 ) 1 : Once again we use induction for the proof. a) k = q = 0 : p0;0 = 0 must hold for all wavelets. b) k + q = l : given that
k+q H 0(!x ; !y ) 0 = @ @! ; for all aliasing frequencies ? !i; and(k + q) l; q xk @!y !=!i
we assume that pk;q = 0; 8p = 1; : : : (M ? 1) ; k + q < l and prove that pk;q = 0 for all k + q = l. The orthogonality condition in the modulation representation (claim 1) gives
H 0(!)H p(!) +
X i
H 0(! + !i )H p(! + !i ) = p;0 ;
where the sum is over all aliasing frequencies. Di erentiating this equation l times, i.e. with (k + q ) = l, we obtain
"
@ k+q H 0(!)H p(!) + X H 0(! + ! )H p(! + ! ) i i @ !kx @ !qy i
Using condition 3,
# = 0: ! =0
! ! st q k X k?s q?t H p(!x ; !y ) X q k @ H ( ! ; ! ) @ x y 0= ; @!xs @!yt ! @!xk?s @!yq?t ! t s s t +
0
+
=0 =0
=0
=0
because of the assumption, only the term of s = t = 0 will contribute to the result, hence we are left with p p @ k+q H p(!x; !y ) M (20) = Mp = 0; p = 1; : : :; M ? 1:
!
@!xk @!yq
k;q
=0
1 ) 3: This can be shown in a very similar manner by starting from the orthogonality of the rows of the modulation matrix, that is
X j
H j (!)H j (! + !k ) = k;0 10
(21)
2
and di erentiating it.
One main di erence between the 2-band case, and M-band wavelets with M 6= 2, is that in the 2-band case the wavelet and scaling functions are fully de ned by each other. The requirement of maximum number of vanishing moments, as can be seen from the third alternative in Theorem 2, implies that the scaling lter resembles an ideal low-pass lter, and in the 2-band case, this implies that the wavelet lter is close to an ideal high-pass lter. While when we have more than one wavelet the requirement of vanishing moments does not give any information about the division of bands between the wavelet lters, and these may overlap extensively. We now deal with the issue of regularity. In the 1D case the relation between the number of vanishing moments of a wavelet function and its regularity has been investigated from several di erent viewpoints[7, 8, 9]. In all three papers cited above, it is concluded that a high number of vanishing moments is a necessary condition for obtaining high regularity, where regularity is measured by the Holder exponent. No such results exist for general N-dimensional wavelets. Furthermore, as suggested in [21, 22], for a given sampling scheme the regularity varies considerably by altering the choice of the sampling matrix and the choice of the coset vectors. The generalization of the 1D derivations is not trivial since they depend, in a crucial way, on the factorability of 1D polynomials. We face this problem, once again, in our attempts to measure the regularity, in terms of the Holder exponent, of a given wavelet function/ lter. A di erent method was suggested by Villemoes, and was adapted to the 2D quincunx case[10]. It works for Holder exponents < 1, and we use it to get a lower bound for the regularity. In the examples considered in the following section, the regularity is assessed by using a simple generalization of the technique suggested in [9]. We de ne the approximation of order-j of the scaling function 0(x) at the points x = D?j n to be : p X (22) h0j (n) = M h0 (n ? Dm)h0j?1 (m); m where h01 (n) = h0(n). The regularity of a scaling lter/function is measured by the rate of convergence of the approximation to the scaling function. Denition 3 The regularity measure of a 2D scaling function/ lter has the value , where N ? 1 < N , if it has at least N vanishing moments and ?j ?N ; 8j ; (23) max max j(l;N ?l)h0j (n)j CM lN n (
2
)
where (k;q) h0j (n) stands for the nite di erence of order (k + q ),
with and
4 (k;q)h0j (n1 ; n2) = [(k?1;q)h0j (n1 ; n2) ? (k?1;q) h0j (n1 ? 1; n2)]=M ?j=2 = [(k;q?1)h0j (n1 ; n2) ? (k;q?1) h0j (n1; n2 ? 1)]=M ?j=2
(1;0)h0j (n1 ; n2) = [h0j (n1; n2 ) ? h0j (n1 ? 1; n2)]=M ?j=2 (0;1)h0j (n1 ; n2) = [h0j (n1; n2 ) ? h0j (n1; n2 ? 1)]=M ?j=2 11
As mentioned above, it is very dicult to give general results regarding the regularity of wavelets in more than 1D. The next theorem, which is a generalization of the 1D version given in [7], and is given for a certain class of sampling matrices, shows, however, that the number of vanishing moments is related to the regularity in more than 1D as well. (Although it is stated for wavelets with compact support, it can be proven for more general wavelets with sucient decay.) Theorem 3 In the case of sampling matrices D, for which there exists an integer n such that Dn = M n I = I , a compactly supported wavelet function (x) that is orthogonal, i.e.
Z
j;n(x) i;m(x) dx = j;i n;m ;
(24)
is in C m+1 only if all its rst (m +1) moments vanish, that is, using the notation of Theorem 2 Z mk;q = xq (x) dx = 0; 8(k + q) = 0; : : :m (25)
Proof: (Based on the 1D proof [23]). First we show that the theorem is correct for m = 0. From the orthogonality, we have for every j j , Z Z 0 = (x) (Dj x ? D j ?j n ) dx = M ?j (D?j + D?j n ) () d: (26) 0
(
0)
0
0
0
If we write j = pn, and take the limit of the above expression as p grows to in nity, we have, Z Z ?pn + D?j n0 ) () d = lim ( ?p + D?j n0 ) () d 0 = plim ( D !1 p!1 =
0
Z
(D?j0 n0) () d = (D?j0 n0 )
Z
0
() d:
(27)
Note that we can change the order of the integral and the limit, since () is in C 1 and is compactly supported. This ends the proof for m = 0, since we can choose the point D?j n0 so that (D?j n0 ) 6= 0. We proceed to prove the Theorem by induction on l. Given that (x) 2 C l+1 , with l m we assume that all moments vanish for (k + q ) = i < l and prove that all the moments with (k + q ) = l, equal zero. We can use Taylor's Theorem. Since () is in C l+1 , denoting = (x; y ) we have Xl 1 Xi i ! @ i (D?j n0) ?pn (k;i?k) ? pn ? j (D ) + Rl (D?pn + D?j n0 ); (D + D n0) = i! k i ? k @x @y i=0 k=0 k (28) where the remanent term - Rl is uniformly bounded and ?j lim Rl ( + D n0) = 0 : 0
0
0
0
0
0
!0
jD?pnjl
Hence, equation (26) can be written as ! i ?j Z i Xl 1 X i @ (D n0 ) (D?pn )(k;i?k) () d 0 = @xk @y i?k i=0 i! k=0 k 0
+
Z
Rl(D?pn + D?j n0 ) () d: 0
12
(29)
Using Dn = I , we have
(D?pn )(k;i?k) = ?pi xk y i?k ; and D?pn i = ?pi jji :
According to our assumption, only the term with i = l is non zero. Hence, by dividing both terms of (29) by ?pl , we have l X l 0 = l1! k=0 k
!
@ l (D?j n0 ) Z xk yl?k () d + Z Rl ( ?pi + D?j n0 ) jjl () d: (30) @xk @y l?k j ?pi jl 0
0
Once more we take the limit of the above equation as p ! 1, and change the order of the integration and limit, to obtain
!
l X l @ l (D?j n0) Z 0 = l1! @xk @y l?k k=0 k ! l ?j Z l X l 1 @ (D n0) = l! k=0 k @xk @y l?k
Z
0
xk yl?k () d +
0
xk yl?k () d + 0:
Rl( lim p!1
?pi + D?j0 n0 ) jjl () d l ? p j j
(31)
In order to complete the proof we still have to show that each term in the above sum is zero. Since (31) is true for every point D?j n0, i.e. for every point in a set dense in IR 2, l this can be done by showing that the l + 1 functions @x@ k @y(x;yl?)k are independent. if 0
0
Xl @ l (x; y) a ;
k=0
k @xk @y l?k
we can Fourier transform this equation and integrate by parts, using the fact that all functions vanish at in nity, we have 0 0
Xl Z 1 @ l (x; y) ?i x!x e a
k=0 Xl
k=0
k
?1 @xk @y l?k
ak !xk !yl?k
Z1
?1
(
y!y ) dx dy ;
+
(x; y )e?i(x!x+y!y ) dx dy:
can not be the zero function, hence, we are guaranteed that there exist no neighborhood at which the Fourier transform of is zero, since is of compact support. Thus the 2D polynomial must be the zero function, i.e. ak = 0 8k, thus proving that the l + 1 functions are independent. 2 It goes without saying that Theorem 3 can be applied to each of the M ? 1 wavelet functions - p(x). Note also, that the proof depends on the orthogonality of the functions j;n(x), and does not hold for wavelets that are only orthogonal in the lter-bank sense, i.e. that are only lossless. The 1D proofs given in [9, 8], use a completely di erent starting point and are given also for non-orthogonal wavelets, but unfortunately have no simple 2D generalization. The speci c examples given in [22, 21] indicate that choosing the sampling matrix so that Dn = I , results indeed in more regular wavelets. 13
4 Designing wavelets with vanishing moments We now turn to the actual design of orthogonal factorable wavelets. A Rank M wavelet matrix is constructed from, paraunitary building blocks of degree one. With the wavelet matrix at hand, the wavelet lters are reconstructed from their polyphase components; (H 0(z1; z2); H 1(z1 ; z2); : : :; H M ?1(z1; z2)) = (z?k ; z?k ; : : :; z?kM ? )HI (zD )H(M ) ; (32) where ki is the shift vector corresponding to the i'th coset. The general form of the Haar wavelet matrix is ! 1 0 H(M ) = H(cM ) (33) 0 OM ?1 where OM ?1 is an orthogonal (M ? 1) (M ? 1) matrix and H(cM ) is the canonical Mband Haar wavelet matrix [17]. The rst column of the canonical Haar wavelet matrix is p1 (1; 1; : : :; 1)T , hence the scaling lter is given by M 0
1
1
0 1 BB 11 CC H (z ; z ) = (z?k ; z?k ; : : :; z?kM ? )HI (zD ) p1 B B .. CCC : MB @.A 0
1
2
0
1
1
(34)
1 Thus a factorable scaling lter is de ned ! by (m + n) unit vectors, that is by (M ? 1)(m + n) M ? 1 parameters which de ne the Haar wavelet matrix free parameters. The other 2 do not e ect the scaling lter, but only the other wavelet lters. Using any of the three alternatives of Theorem ! ! 2, one can see that for L vanishing L +1 moments we have Nc = (M ? 1) ? 1 constraints, hence in general we need 2 wavelets with at least Nfp = Nc free parameters. If one uses the third alternative of Theorem 2, that is equations (17), it is obvious that only the parameters de ning the scaling lter can be taken into account. The other parameters do not appear in equation (34) and therefore do not appear in the constraint equations (17) either. This decoupling of the Haar wavelet parameters from the vanishing moment constraints is the reason for de ning the polyphase (wavelet) matrix the way we did, instead of the more common way which uses the transpose of the one we use. (Another way of reaching this goal is to use a di erent factorization, in which the characteristic matrix is the rst factor on the left, and not on the right.) These considerations suggest that in order to obtain wavelets with ! L vanishing moments we should L + 1 ? 1. Furthermore, in both cases use wavelets with McMillan degrees (m + n) 2 dealt with, the two-band quincunx and four-band separable sampling, it is obvious from 1D results that each polynomial degree should be larger or equal to L ? 1. Otherwise we could nd 1D projections with L vanishing moments and length less than 2L, contrary to 1D well known results [5]. The rest of the section is devoted to two types of orthogonal 2D wavelets, quincunx and four-band separable sampling wavelets. 14
4.1 Quincunx Wavelets
The quincunx sampling is a 2D, two-band, nonseparable sampling scheme. Detailed analysis of 2D multirate systems with quincunx sampling can be found in Kovacevic's doctoral thesis [3]. The sampling and Haar wavelet matrices we use are (there are many other possibilties):
!
1 1 ; 1 ?1
D=
H = p12 (2)
!
1 1 ; 1 ?1
and choose ki = f(0; 0)T ; (1; 0)T g as the coset shift vectors. With this choice of the sampling matrix we have D2 = 2I , and it has been shown to result in wavelets with higher regularity than with other sampling matrices[3, 22, 21]. In the quincunx case the McMillan degrees (m; n) of the wavelet matrix are equal to the polynomial degrees of the elements of the matrix, so for the rest of the subsection we simply use the term degree. Similar to the two-band 1D case, the wavelet lter is fully determined by the scaling lter. For a scaling lter given by: H 0(z1; z2) = H00(z1z2 ; z1z2?1 ) + z1?1 H10(z1z2 ; z1z2?1 ); the wavelet lter is:
H 1(z1 ; z2) = (z1 z2)?m (z1z2?1 )?n H10 (z1?1z2?1 ; z1?1z2 ) ? z1?1(z1z2)?m(z1z2?1)?n H00(z1?1z2?1; z1?1z2): In the two-band factorable case the vectors Vk in equation (9) are 2D unit vectors and thus can be parameterized by one parameter, V (k ) = (cos(k ); sin(k ))T , where k 2 [0; 2 ): The general form for the vanishing moment constraints on the angles k that de ne the wavelet system is developed in Appendix A. We use the third alternative in Theorem 2, i.e. , we require the frequency response of the scaling lter to have a zero ! of order L at the L + 1 ? 1 constraints aliasing point (z1 = ?1; z2 = ?1). We get the following Nc = 2
p+l 0 0 = @ Hp (z1l; z2) @z @z X 1 2 m;nz =z =?1 = ?2?k;q (p; l)Sk;q + p?m;n k;q (p ? 1; l)(Sk;q + Ck;q ) 1
k q even
X
( + )
+
2
k q odd
( + )
m;n 2?m;n k;q (p; l)Ck;q + p?k;q (p ? 1; l)(Sk;q ? Ck;q );
(35)
with (p + l) = 1; 2; : : :; (L ? 1). The coecients ?m;n k;q (p; l) were computed and tabulated, +n . As mentioned earlier, and Sk;q and Ck;q are trigonometric functions of the angles fgmk=1 in order to achieve L vanishing moments one should use in general wavelets of degrees (m + n) Nc . Since the wavelets with vanishing moments are obtained by solving a set of non-linear equations, in general we will have more than one solution with minimal degrees, for a given number of vanishing moments. This is true for the 1D case as well, but contrary 15
to the 1D wavelets with vanishing moments dealt by Daubechies, here di erent solutions may have completely di erent spectral characteristics. We would like to mention the special case of solutions that are actually 1D lters.
Lemma 2 A degree L?1 1D wavelet lter with L vanishing moments is likewise a factorable quincunx wavelet lter with degrees (L ? 1; L ? 1) and L vanishing moments. This can easily be demonstrated by factorizing each degree-1 factor of the 1D polyphase matrix as
I + (z12 ? 1)ViViT = I + (z1z2 ? 1)ViViT I + (z1z2?1 ? 1)ViViT :
These 1D quincunx wavelets were suggested by Cohen and Daubechies[22]. Although the wavelet lters are 1D, due to the non-separable sub-sampling the resulting wavelet functions are 2D non-separable functions. Thus, the above arguments concerning the minimal degrees needed in order to obtain L vanishing moments are true in general, but for the special case of 1D lters we can satisfy the constraints, using wavelets of degrees (L ? 1; L ? 1) even though 2(L ? 1) < Nc for all L > 2. We now give a few examples, but rst note that:
The vectors Vi appear in H(z ; z ) (see equation (9)) only in the combination ViViT , 1
hence, since
2
V (i )V (i )T = V ( + i )V ( + i)T we can limit the search to i 2 [? 2 ; 2 ). for every solution fig there exists another solution f( ? i)g, that will give rise to a scaling lter (function) which is the re ection of the original lter (function) through its center, i.e. , ~h0 (n1 ; n2) = h0 (N1 ? n1 ; N2 ? n2 ).
Given a solution h (n ; n ) with degrees (m; n) there exists another solution with degrees (n; m) in which the variable z is interchanged with z in the polyphase matrix. This will cause the transformation hi (n ; n ) ! hi (n ; ?n ). Hence the following examples are restricted to degrees (m; n) where m n. Example 1.1: The smallest quincunx wavelet lter is the Haar wavelet, its degrees are (m; n) = (0; 0), thus the corresponding wavelet matrix is the Haar wavelet matrix H . 0
1
2
1
2
1
2
1
2
(2)
The lters are exactly the 1D Haar lters, i.e. h0(n1 ; n2) := p1 (1 ; 1) and h1(n1; n2) := p1 (1 ; ?1); 2 2 and only the rst moment vanishes. The absolute value of the frequency response of the scaling lter, and the scaling function in this case are plotted in Fig. 1.a and 1.b respectively. Example 1.2: We next deal with wavelets with L = 2 vanishing moments, in this case the minimal degrees are (1; 1). The two constraint equations are: 1.
0 = 1 ? cos(21 ) ? cos(22) 16
2. 0 = cos(21 ) ? cos(22): We have two di erent solutions (and their re ections), 1. 21 = 22 = 3 : which in accordance with Lemma 2 is the 1D, 4-tap lter of Daubechies with two vanishing moments.
1 0 0 0 0 0 B p p p p C h (n ; n ) := p1 B 1+ 3 3+ 3 3? 3 1? 3 C A @ 4 2 0
1
2
0
0
0
0
2. 21 = ?22 = 3 : which produces the scaling lter
1 0 p p 3 3 + 3 0 0 3 + 3 B p p p p C: h (n ; n ) := p1 B ? 1 ? 3 3 + p3 3 ? p3 ?1 + 3 C A @ 8 2 0 0 3? 3 3?3 3 0
1
2
The frequency response of the scaling lter and the scaling function for the rst solution are shown in Fig 1.c and 1.d, and the corresponding functions for the second solution are depicted in Fig. 1.e and 1.f. For L > 2 the constraint equations can not be solved analytically and a numerical program given in Numerical Recipes [24] was used for the solution of the non-linear equations. ! m + n In addition to the (m + n) free parameters de ning the wavelet matrix, there are
n
alternatives in which the factors of degree-1 can be arranged in equation (9). Each alternative results in a di erent set of constraints, thus introducing another discrete parameter into the problem. It turns out that a solution does not exist for all possible alternatives. For example, analyzing the equations for wavelet lters with degrees (3; 2) we found that one can not obtain three vanishing moments in the two cases in which all the factors of degree-1 in one variable precede the factors in the other variable. Example 1.3: The minimal degrees for obtaining three vanishing moments (not including 1D lters) are (3; 2). The functions plotted in Fig. 2.a and 2.b correspond to a solution obtained by choosing the combination (z1; z1; z2; z1; z2) and the ve parameters as
f2gi = f?0:896794; ?2:934522; 0:782106; 0:546171; 1:781789g: 5 =1
In Figures 2.c and 2.d we show another solution, this solution was obtained using the order (z1 ; z1; z2; z2; z1), and the parameters
f2gi = f0:808218; 1:258868; 0:782106; 1:781789; ?2:091721g: 5 =1
Although both solutions have three vanishing moments there is a considerable di erence between them. The rst solutions is much closer to an ideal quincunx low-pass lter, with 75% of its energy concentrated in the low-pass band, as opposed to the second solution with 60%. The di erence becomes also apparent by looking at the scaling functions. The 17
regularity measure, introduced in the previous section, gave the value = 0:6 0:1 for the rst solution, and = 0:05 0:1 for the second. We calculated the lower bound of for the second example. Using the method suggested by Villemoes[10], we were not able to show that the function is even continuous ( > 0), when applying the method for resolution levels j 14. The 1D solution, corresponding to the Daubechies wavelet with three vanishing moments, obtained the highest value of regularity with = 1:0 0:1, which agrees with the 1D bound for the Holder exponent of = 1:0878. For the sake of completeness we searched the entire space of wavelet matrices with degrees (3; 2) and three vanishing moments, we did so using the complete charicterazation of 2D PUP matrices given in [12]. To obtain compactly supported wavelets with three vanishing moments, one must solve 15 equations for 15 variables. These equations were solved numerically, using the same procedure mentioned above. All the solutions found coincide with solutions that were obtained using factorable wavelets, but it is feasible that other solutions, that were not found numerically, do exist. Finally we give two examples of quincunx scaling lters and functions with four vanishing moments. The Frequency response of the scaling lter and the scaling function are plotted in Fig. 3.a and 3.b for one example and in Fig. 3.c and 3.d for the other. The rst example has a good frequency response with 80% of the energy in the low-pass band, but is very spiky, = 0:5 0:1. While the second example is more regular, with = 0:9 0:1, but with a slightly inferior frequency response.
4.2 Four-band, separable sampling wavelets
This sampling scheme is the most commonly used in all applications of the 2D DWT , since it is the simplest ! generalization of the two-band 1D sampling. The sampling matrix is 2 0 D= and the four-band Haar wavelet matrix has the general form; 0 2
0 1 p 1 ? 3 0 0 q B CC p ! B 3 ? 2 0 1 1 = B C 1 0 1 q B C p ; p C H = 2B 2 0 O 1 1 = 3 ? B C @ p q p A 2 3
(4)
1 1= 3
2 3 2 3
3
(36)
2
O3 being a rank-3 orthogonal matrix. Thus, we have a three parameter space of Haar wavelet matrices. We chose the four coset shift vectors to be: ki = f(0; 0)T ; (1; 0)T ; (0; 1)T ; (1; 1)T g, hence, the three points of repeated spectra are: ! i = f(; 0)T ; (0; )T ; (; )T g: With this choice the expression for a factorable scaling lter of McMillan degree (m; n) is according to (34):
0 1 1C B m n B C Y H (z ; z ) = 21 (1; z ; z ; z z ) [I + (zk? ? 1)Vk VkT ] B BB 11 CCC : k @ A 0
+
1
2
1
2
1 2
=1
18
2
1
(37)
Once again, in order to obtain L vanishing moments we impose Nc = 3
!
L+1 ?1
!
2 constraints on the 3(m + n) free parameters de ning the scaling lter. Just as we did in the quincunx case, we seek solutions with the minimal ! possible degrees which in general are L + 1 ? 1. Once we have designed the scaling given by Nfp = Nc , that is (m + n) = 2 function by xing the 3(m + n) free parameters we still have three free parameters de ning the Haar wavelet matrix. Note that using 1D wavelets with L vanishing moments to design separable 2D wavelets will result in wavelets with McMillan degrees ! (2L ? 2; 2L ? 2) (see section 2.2). Since in L + 1 ? 1 for all L 6= 6, in the following examples we this case (m + n) = 4(L ? 1) 6= 2 will not obtain separable lters as private solutions. Example 2.1: The Harr wavelet matrix gives rise to the smallest wavelet lters. Each lter has four non-zero coecients, and the corresponding function is compactly supported on [0; 1]2. In this case only the rst moment vanishes and the resulting functions are discontinuous and piecewise constant. Example 2.2: In order to obtain wavelets with L = 2 vanishing moments we must use wavelets with McMillan degrees (1,1), thus the wavelet matrix is given by,
H(z ; z ) = [I + (z? ? 1)V V T ][I + (z? ? 1)V V T ]H : 1
2
1
1
1 1
2
1
2 2
(4)
We have six bilinear constraints, (two derivatives at three aliasing frequencies) on eight unknowns, (two vectors of length 4) de ned by six free parameters. The two unit vectors solving this set of equations are, p p p p V1 = p1 ( 3; 1; 3; 1)T ; V2 = p1 ( 3; 3; 1; 1)T : 2 2 2 2 As mentioned previously, the three parameters de ning the Haar wavelet matrix do not in uence the number of vanishing moments but rather control, in some sense, the orientation of the three wavelet lters/functions. In the following example we used the separable Haar wavelet matrix 1 0 1 1 1 1 C BB 1 ?1 1 ?1 C 1 C B (4) H = 2B B@ 1 1 ?1 ?1 CCA ; 1 ?1 ?1 1 this is one way to guarantee that each wavelet lter has a maxima at one aliasing point and a zero at all others. Once again, for L > 2 only numerical solutions can be obtained. In Fig. 4 an example of wavelet lters and functions with four vanishing moments are depicted. Their McMillan degrees are (5; 4), the lters are of length (12; 10), and the regularity of the scaling function was measured to be = 1:25 0:1. 19
5 Discussion and Conclusions Several methods have been suggested for constructing 2D, orthogonal wavelet lters and functions [1, 22]. These methods are based on some sort of transformation from 1D wavelets. Hence, the resulting wavelets in each case, constitute only a small restricted subset of all possible orthogonal wavelets. Here we use a parameterization of all factorable orthogonal wavelets, which is a much more general subset. The subset of factorable wavelets was shown in [12] to include almost all possible orthogonal wavelets, and contrary to the complete set of orthogonal wavelet, has a simple parameterization. We have shown that this subset includes the previously suggested 2d orthogonal wavelets, such as, separable wavelets and the special (actually 1D) quincunx wavelets [22]. In order to obtain wavelets with several vanishing moments, one has, in general, to solve numerically a set of non-linear equations. Even when focusing our attention only on wavelets with maximal number of vanishing moments for given McMillan degrees, considerable freedom is still left, since the set of non-linear equations have many solutions. Furthermore, usually, satisfying the conditions of having a large number of vanishing moments is not enough, as is illustrated in Fig. 2 and Fig. 3, Extra care should be taken in choosing the best solution for a given 2D application. Theorem 4 and our examples show that, a sucient number of vanishing moments is a necessary, but by no means sucient, condition for the regularity of the wavelet function. The regularity of the original Daubechies orthogonal wavelets (minimal phase wavelets) was found to grow at least linearly with the number of vanishing moments, yet in 2D we were not able to select any special subset of solutions and to give a lower bound for their regularity. If a high degree of regularity is required for given McMillan degrees, one should probably sacri ce some of the vanishing moments and use the extra degrees of freedom to optimize the regularity. For a given order of moments, one can relax only part of the partial derivative constraints (17), and only at some of the aliasing points. Our treatment of the regularity issue is quite preliminary and mainly serves as a motivation for investigating 2D wavelets with vanishing moments. Further work is required in order to establish the exact relation between the number of vanishing moments and regularity, and its dependency on the choice of sampling matrix and coset vectors. As in many other areas of mathematics, most of the diculties encountered when trying to generalize 1D theory to higher dimensions arise already at the transition from 1D to 2D. Once these diculties have been overcome, the design of N-dimensional wavelets will be a simple generalization of the 2D method.
Appendix A In this appendix we will develop the constraint equations for vanishing moments in the case of factorable quincunx wavelets. Each factor of degree-1 is de ned by one parameter i
Wi = (I + (zi? ? 1)ViViT ); 1
20
were zi = z1 or z2 and Vi is a unit vector of length 2. If we use Vi = (cos(i ); sin(i ))T ; we get !! cos(2 ) sin(2 ) 1 i i (A:1) Wi(i) = 2 (zi?1 + 1)I + (zi?1 ? 1) sin(2 ) ? cos(2 ) : i i Next we de ne two matrices, the 2D rotation matrix- R() and the rotation-plus-re ection matrix T () , ! ! cos(2) ? sin(2) cos(2) sin(2) ; T () = ; R() = sin(2) cos(2) sin(2) ? cos(2) and recall the relations: T ( )T ( ) = R( ? ); R( )R( ) = R( + ); R( )T ( ) = T ( + ); T ( )R( ) = T ( ? ): (A:2) From this we see that the 2D polyphase matrix of McMillan degrees (m; n) can be written as : mY +n m X n X m;n(z ; z )T (fg); H(z1; z2) = W(i) = (A:3)
k;q 1 2 k;q k=0 q=0
i=1
where we de ne 4 m;n(z ; z ) =
k;q 1 2
1 m
n
+
(z1?1 + 1)m?k (z1?1 ? 1)k (z2?1 + 1)n?q (z2?1 ? 1)q ; 2 and Tk;q (fg) is the sum of all possible combinations of an ordered product of k T -matrices of z1 factors and q T -matrices of z2 factors. With use of the relations (A.2) one can see that for a given k and q , Tk;q (fg) will include terms of the following type
8 > > > 4< Tk;q (fg) = > > :
! ?Sk;q =4 P R( ? + ? ? ); if (k+q) is even j j j jk q Ck;q ! Sk;q =4 P T ( ? + ? + ); if (k+q) is odd j j j jk q ?Ck;q For example with degrees (m; n) = (2; 2) ,where the rst two factors are in z and the last in z , we have 16H(z ; z ) = (z ? + 1) (z ? + 1) + (z ? ? 1)(z ? + 1) [T ( ) + T ( )] + (z ? + 1) (z ? ? 1)[T ( ) + T ( )] + (z ? ? 1)(z ? ? 1)[R( ? ) + R( ? ) + R( ? ) + R( ? )] + (z ? ? 1) (z ? + 1) [R( ? )] + (z ? + 1) (z ? ? 1) [R( ? )] + (z ? ? 1) (z ? ? 1)[T ( ? + ) + T ( ? + )] + (z ? ? 1)(z ? ? 1) [T ( ? + ) + T ( ? + )] + (z ? ? 1) (z ? ? 1) [R( ? + ? )]: Ck;q Sk;q Ck;q Sk;q
1
2
3
+
1
2
3
+
1
2
1
2
1
1
1
1
1
1
1
1
1
1
2
2
1
2
2
2
2
1
2
1
2
1
2
2
1
2 2 2
2
2
2
2 1
2 1 2
1
2
2
2
1
2
3
4
1
1
2
1
2
2
2
1
2
3
1
1
2
3
4
4
2
3
1
2
3
1
2
4
1
3
4
2
3
4
1
2
21
3
4
2
4
With an explicit expression for H(z1; z2) we can obtain the vanishing moments constraints as constraints on the parameters . First we construct the scaling lter from the polyphase matrix ! 1 1 ; (A:4) H 0(z1 ; z2) = (1; z1)H(z1z2 ; z1z2?1 ) p 2 1 and then we write the partial derivatives of the scaling function
p @ p lH (z ; z ) 2 @z p @z l z +
0
1
m X n h X k=0 q=0
1
2
2
z ?1
=
1= 2=
i
m;n ?m;n k;q (p; l)(Ck;q(1 1) ? Sk;q (1 1)) + p?k;q (p ? 1; l)(Sk;q Ck;q ) ;
where the upper sign should be used when (k + q ) is even and the lower sign otherwise, and we used the notation m;n ?1 4 @ p+l k;q (z1 z2 ; z1z2 ) m;n : ?k;q (p; l) = p l
@z1 @z2
z1 =z2 =?1
Finally we arrive at the constraints for L vanishing moment. ! By equating all partial derivaL +1 tives (p + l) = 1; 2; : : :; (L ? 1) to zero we get ? 1 constraints 2
i X h m;n ?2?k;q (p; l)Sk;q + p?m;n ( p ? 1 ; l )( S + C ) k;q k;q k;q k q even i X h m;n + 2?k;q (p; l)Ck;q + p?m;n k;q (p ? 1; l)(Sk;q ? Ck;q ) :
0 =
( + )
k q odd
(A.5)
( + )
Acknowledgments
We would like to thank N. Cohen, M. Zibulski of the Department of Electrical Engineering, and Y. Binyamini of the Department of Mathematics, Technion, for their help. We also wish to thank L. Villemoes for his help in applying his methods to our examples.
References [1] S. Mallat, \A theory of multiresolution signal decomposition: The wavelet representation," IEEE Trans. PAMI, vol. 11, pp. 674{693, 1989. [2] T. Chen and P. P. Vaidyanthan, \Mulidimensional multirate lters and lter banks derived from one-dimensional lters," IEEE trans. on SP, vol. 41, May 1993. [3] J. Kovacevic, Filter banks and Wavelets: Extentions and Applications. PhD thesis, Department of Electrical Engineering, Columbia University, 1991. 22
[4] C. Lien and C. Huang, \The design of 2-d nonseparable directional oerfect reconstruction lter banks," Multidimensional Systems and Signal Processing, vol. 5, p. 289, 1994. [5] I. Daubechies, \Orthonormal bases of compactly supported wavelets," Comm. Pure and Appl. Math, vol. 41, pp. 909{996, 1988. [6] O. Rioul, \On the choice of wavelet lters for still image compression," in ICASSP, vol. V, p. 550, 1993. [7] I. Daubechies, Ten lectures on wavelets. CBMS/NSF series in applied math., SIAM, 1992. [8] I. Daubechies and J. C. Lagarias, \Two-scale di erence equations II. Local regularity, in nite products of matrices and fractals," SIAM J. Math. Anal., vol. 23, no. 4, pp. 1031{1079, 1992. [9] O. Rioul, \Simple regularity criteria for subdivision schemes," SIAM J. Math. Anal., vol. 23, no. 6, pp. 1544{1576, 1992. [10] L. F. Villemoes, \Continuity of nonseparable quincunx wavelets," Appl. Comput. Harmon. Anal., vol. 1, pp. 180{187, 1994. [11] G. Karlsson and M. vetterli, \Theory of two dimensional multirate lter banks," IEEE trans. on ASSP, vol. 38, June 1990. [12] S. Venkataraman and B. C. Levy, \State space representations of two dimensional FIR lossless transfer matrices," IEEE Trans. Circuits and Syst., vol. 41, p. 117, Feb. 1994. [13] S. Venkataraman and B. C. Levy, \A comparison of design methods for 2-D FIR orthogonal perfect reconstruction lter banks." submitted to IEEE Trans. Circuits and Syst., 1994. [14] J. Kovacevic and M. Vetterli, \Nonseparable two and three-dimensional wavelets," IEEE Trans. Signal Process., vol. 43, May 1995. [15] E. Viscito and J. P. Allebach, \The analysis and design of multidimensional FIR perfect reconstruction lter banks for arbitrary sampling lattices," IEEE Trans. Circuits and Syst., vol. 38, pp. 29{41, Jan. 1991. [16] W. M. Lawton and H. L. Resniko , \Multidimensional wavelet bases," tech. rep., Aware, Aware Inc. Cambridge MA, 1991. [17] P. N. Heller, H. L. Resniko , and R. O. Wells, \Wavelet matrices and the representation of discrete functions," in Wavelets: a tutorial in theory and applications (C. K. Chui, ed.), Academic press, 1992. [18] P. P. Vaidyanthan, T. Q. Nguyen, Z. Doganata, and T. Saramaki, \Improved technique for design of perfect reconstruction FIR QMF banks with lossless polyphas matrices," IEEE Trans. on ASSP, vol. 37, July 1989. 23
[19] D. Stanhill and Y. Y. Zeevi, \Two-dimensional orthogonal wavelets, with linear phase," Tech. Rep. CC PUB#111, Dep. of Electrical Eng. , Technion, Haifa, ISRAEL, 1995. [20] P. Ste en et al., \Theory of regular m-band wavelet bases," Tech. Rep. TR-93-02, Computational Math. Lab., Rice University, Houston, TX, 1993. [21] K. Grochenig and W. R. Madych, \Multiresolution analysis, Haar bases, and selfsimilar tilings of Rn ," IEEE Trans. Inform Theory, vol. 38, March 1992. [22] A. Cohen and I. Daubechies, \Non-separable bidimensional wavelet bases," Rev. Mat. Iberoamericana, vol. 9, no. 1, pp. 51{137, 1993. [23] G. G. Walter, Wavelets and other orthogonal systems with applications. CRC Press, 1994. [24] W. H. Press et al., Numerical Recipes in C. Cambridge University Press, 1992.
FIGURE CAPTIONS Fig. 1: Quincunx scaling lters and functions. The frequency response of the scaling
lter, and the scaling function of: (a) and (b) the Haar quincunx wavelet (Example 1.1), (c) and (d) wavelet with two vanishing moments, 1D solution (Example 1.2), (e) and (f) wavelet with two vanishing moments, second solution (Example 1.2).
Fig. 2: Quincunx scaling lters and functions with three vanishing moments. The
frequency response of the scaling lter, and the scaling function of the rst ((a) and (b)) and second ((c) and (d)) solutions of Example 1.3. Both scaling functions show the eighth order approximation, using the sub-division scheme.
Fig. 3: Quincunx scaling lters and functions with four vanishing moments. The frequency response of the scaling lter, and the scaling function of the two examples mentioned in the text.
Fig. 4: Four-band, separable sampling wavelet with four vanishing moments. Figures (a)-(d) show the magnitude response of the four wavelet lters, and (e)-(h) show the wavelet functions.
24
