Wavelet transforms for discrete-time periodic signals - Semantic Scholar
SIGNAL PROCESSING, VOL. 42, NO. 2, PP. 167{180, 1995
Wavelet transforms for discrete-time periodic signals
John A. Gubner and Wei-Bin Chang
Department of Electrical and Computer Engineering, University of Wisconsin, Madison, WI 53706, USA
Abstract. Wavelet transforms for discrete-time periodic signals are developed. In this nite-dimensional context, key ideas from the continuous-time papers of Daubechies and of Cohen, Daubechies, and Feauveau are isolated to give a concise, rigorous derivation of the discrete-time periodic analogs of orthonormal and symmetric biorthogonal bases of compactly supported wavelets. These discrete-time periodic wavelets are expressed in terms of circular FIR lters, and thus lead to fast wavelet transforms whose complexity is order N . Keywords: Orthonormal wavelets; Biorthogonal wavelets; Symmetric biorthogonal wavelets; Pyramid algorithms; Fast wavelet transforms
1. Introduction
tion 2 we review the properties of abstract pyramid-type algorithms [4], which are the key to fast wavelet transforms. Section 3 specializes these abstract pyramid algorithms to discrete-time periodic signals, de nes the discrete wavelet transform, and discusses fast wavelet transforms. Section 4 derives the discrete-time periodic analogs of Daubechies' orthonormal bases of compactly supported wavelets [9]. Section 5 derives the discrete-time periodic analogs of the symmetric, biorthogonal bases of compactly supported wavelets in [7]. 1.1. Sinusoids and wavelets Let X0 denote the space of discrete-time signals x0 (n) with period N, where N is a power of 2. Before de ning wavelets, recall that from the theory of discrete Fourier transforms, every signal x0 2 X0 can be expressed in terms of N complex sinusoids as
Although there is a large and growing literature on wavelets, there is, to our knowledge, no reference that focuses on the development of wavelets for discrete-time periodic signals, and that includes a derivation of Daubechies' FIR lter coecients in this nite-dimensional context. This paper is an attempt to ll this gap, which is making its presence felt as researchers use discrete-time periodic wavelets in practical applications ranging from image compression [17, Section 13.10] to numerical analysis [1]{[3], [8]. Because discrete-time periodic signals constitute a nitedimensional vector space, our presentation is based primarily on concepts from linear algebra and discrete Fourier transforms, though some details in Appendix A require results about the factorization of polynomials. Hence, we are able to give a concise, rigorous exposition of the key ideas in [9] and [7] that pertain to the discrete-time periodic analogs of compactly supported orthonormal and symmetric biorthogonal wavelets. Our presentation is aimed at the non-wavelet-specialist who is not familiar with the more technically demanding continuous-time papers [9] and [7]. The relationship between the continuous wavelet transform and the discrete wavelet transform for discrete-time aperiodic signals has been studied by Shensa [21]. Multiresolution analysis for continuous-time signals was developed by Mallat [15], and multiresolution analysis for discrete-time aperiodic signals was developed by Rioul [18]. Wavelet transforms for nite-length sequences have been developed by Caire et al. [5], and some of their equations can be found in our paper. However, [5] focuses on signals taking values in nite elds and does not consider the consistency of Daubechies' constraint (37). Also, [5] does not consider biorthogonal wavelets for discrete-time periodic signals as we do in Section 5. The connection between wavelets and lter banks has been well documented by Vetterli and Herley [26]. Recent expository treatments of wavelets include [6], [10], [12], [13], [20], [23], and [24]. The paper is organized as follows. In the remainder of this section, we introduce discrete-time periodic wavelets and compare them to discrete-time sinusoids. In Sec-
x0(n) = where
j2kn=N
NX ?1 k=0
hx0; Ek iEk (n); NX ?1
Ek (n) = e p ; hx0 ; Eki = x0 (n)Ek (n); N n=0 and the overbar denotes the complex conjugate. The rea?1 son the decomposition works is that the signals fEk gNk=0 form an orthonormal basis for X0 . While direct computation of the Fourier coecients hx0 ; Eki for k = 0; : : :; N ? 1 would require N 2 multiplications, the fast Fourier transform can compute these coecients with order N log2 N multiplications. In this paper we consider orthonormal \wavelet" bases for X0 . To create a wavelet basis of dimension N, we rst choose the number of scales or resolution levels, denoted by J, where N=2J 2. The basis is then partitioned according to resolution level as follows. At level j, j = 1; : : :; J, there are N=2j wavelet signals denoted by j;0; : : :; j;N=2j ?1 . As this accounts for only N ? N=2J signals, the basis is completed by adding N=2J scaling functions, 'J;0 ; : : :; 'J;N=2J ?1 . Since these signals are chosen 4
1
4
to form an orthonormal basis, every x0 2 X0 has the ex- uct on U V by taking pansion u ; u~ = hu; u~i + (v; v~): v v~ 2j ?1 J N=X X hx0 ; jii ji (n) x0(n) = Now, if X is a third inner-product space, and if G: X ! U j =1 i=0 (1) and H: X ! V are linear transformations, we de ne the N=X 2J ?1 mapping T: X ! U V by hx0; 'Ji i'Ji (n): + G x = Gx : Tx = H (3) i=0 Hx The numbers hx0; jii and hx0 ; 'Jii are called wavelet co- Next, if the adjoints G : U ! X and H : V ! X exist, ecients, and the fast wavelet transform is an algorithm then the adjoint T : U V ! X also exists and is given by for computing them with order N multiplications. Ex u u amples of wavelets for the case N = 64 and J = 4 are T v = [G H ] v = Gu + H v: shown in Figs. 1{5. It will be shown later that for each j, ji(n) = j 0 (n ? i2j ) for i = 0; : : :; N=2j ? 1. In Our interest is in isometries; i.e., operators that preserve other words, the ji are just circular translates of j 0 inner products, or equivalently, operators for which T T = by multiples of 2j . Similarly, 'Ji (n) = 'J 0(n ? i2J ) for I. For such operators, x can obviously be recovered from i = 0; : : :; N=2J ? 1. We also point out that a quick glance Tx. at the gures might lead one to wonder if in analogy with Remark. If T T = I and T is onto, then TT = I as continuous time, well [11, p. 185]. Such an operator is said to be unitary. Of ( p course, if dim X = dim U V < 1, then T T = I always 2 j +1;2(2n); j +1;2(2n) 6= 0; ? (2) implies TT = I, since in this case T is 1-to-1 if and only j;2(n) = 0; otherwise: if T is onto [14, p. 81]. However, a close inspection of Figs. 1 and 2 reveals that For future reference, note that since T has the form in this cannot be true since 2;2 has 10 nonzero components, (3), T T = I is equivalent to while 1;2 has only 4 nonzero components. To understand G G + H H = I; (4) more fully what is happening, it is necessary to consider the limiting continuous-time case [9], [10]. and TT = I is equivalent to the following three equations: We now compare and contrast the sinusoidal basis funcGG = I; (5) tions Ek with thepwavelets ji and 'Ji . The constant func HH = I; (6) tion, E0(n) 1= N, exhibits no variation at all, while the signals 'Ji (n) exhibit the least variation of all the wavelet GH = 0: (7) signals. At the other extreme, the signal EN ?1(n) exhibits To summarize, if we \decompose" x with the most variation of all the sinusoids Ek (n), while the signals 1;i(n) exhibit the most variation of all the wavelet u = Gx; v = Hx; signals. One important di erence between the wavelets in and, if (4) holds, then we can reconstruct x from u and v Figs. 1{5 and the sinusoids Ek is that the ji(n) are zero with for most values of n. For example, in Fig. 1, we see that x = Gu + H v: 1;2 (n) is nonzero only for n = 2; 3; 4; 5. Thus, looking at wavelet coecients can help localize high-variation features 2.2. The general algorithms in time rather precisely. Slower-varying signals cannot be Let X0 ; X1; : : :; XJ and V1 ; : : :; VJ be two sequences localized as well. of inner-product spaces. Denote the inner product on Xj by h; ij and the inner product on Vj by (; )j . For 2. Pyramid algorithms j = 0; : : :; J ? 1, let Gj : Xj ! Xj +1 and let Hj : Xj ! Vj +1 The basic building block for pyramid-type algorithms be linear operators (for example, see the beginning of Secis a linear transformation from an input space X into a tion 3). Consider the following procedure. product space, say U V , as discussed in Section 2.1. In Section 2.2 we iterate a sequence of these transformations Decomposition Algorithm to obtain a Decomposition Algorithm and a Reconstruction Let x0 2 X0 be given. Algorithm. for j = 0 to J ? 1 xj +1 = Gj xj 2.1. Linear algebra preliminaries vj +1 = Hj xj Let U and V be inner-product spaces with respective next j inner products h; i and (; ). This induces an inner prod- end 4
2
If (4) holds for each pair Gj and Hj , then x0 can be the N-length discrete Fourier transform (DFT) of a se?1 x(n)e ?j2nu=N , recovered from v1 ; : : :; vJ ; xJ by the following procedure. quence x 2 Cj N is de ned by x^(u) = PNn=0 for all integers u. Reconstruction Algorithm De nition (Discrete wavelet transform). For j = Let v1 2 V1 ; : : :; vJ 2 VJ ; xJ 2 XJ be given. 0; : : :; J ? 1, if the linear operators Gj : Xj ! Xj +1 and for j = J ? 1 to 0 H : Xj ! Vj +1 have the form j xj = Gj xj +1 + Hj vj +1 next j N=X 2j ?1 end gj (n ? 2m)x(n); (11) (Gj x)(m) = n=0 If we let N=X 2j ?1 hj (n ? 2m)x(n); (12) (Hj x)(m) = Y = V1 V J XJ ; 4
n=0
and if v1; : : :; vJ , and xJ are generated by the Decomposi- where g , h , and x are N=2j -periodic sequences, and where j j tion Algorithm, then setting gj and hj are such that (4){(7) hold for Gj and Hj , then Dx0 = (v1; : : :; vJ ; xJ ) (8) Dx0 de ned by using these operators in the Decomposition Algorithm is called a discrete wavelet transform (DWT) of de nes a mapping D: X0 ! Y . We equip Y with the in- x0. duced inner product, denoted by (; ). The induced norm We rst note that (11) and (12) are just subsampled is denoted by k k. Note that we are now using row-vector circular convolutions of x(n) with gj (?n) and hj (?n), renotation for product-space vectors. spectively. Next, it can be seen from Section 4.1 and from With the preceding setup, the following facts are easily (57) in Appendix A, that if Gj and Hj satisfy (4){(7), then veri ed. gj and hj are discrete-time periodic analogs of the conju(1) Assuming only that Gj and Hj exist, the vector x~0 gate quadrature lters in [22]. In Section 4 we construct generated by applying the Reconstruction Algorithm gj and hj under two additional constraints. The rst is to an arbitrary vector y~ = (~v1 ; : : :; v~J ; x~J ) 2 Y , is that gj (n) and hj (n) be circular FIR lters, i.e., they are equal to D y~; i.e., the Reconstruction Algorithm im- zero for most values of n in one period. The second is that plements D . Daubechies' constraint given below in (37) hold, forcing gj (2) If (4) holds for each pair Gj and Hj , then D D = I; to be a low-pass lter with maximally at stopband. It i.e., D is an isometryPand preserves norms, and by (8), then follows that on account of (7), hj will be a high-pass kx0 k20 = kDx0 k2 = Jj=1 kvj k2j + kxJ k2J . lter. When gj is a low-pass lter and hj is a high-pass (3) If (5){(7) hold for each pair Gj and Hj , then DD = I. lter, we see that the Decomposition Algorithm rst takes If (4) also holds, then D is unitary. the N-periodic input x0 and computes x1 and v1 , which (4) If D is unitary, then we can write X0 as a direct sum are both N=2-periodic. Clearly, x1 is the result of low-pass of orthogonal subspaces, X0 = W1 WJ BJ , ltering x0 and then subsampling. Similarly, v1 is the rewhere sult of high-pass ltering x0 and then subsampling. The next time through the loop, x1 is low-pass and high-pass Wj = fD y : y = (0; : : :; 0; vj ; 0; : : :; 0); vj 2 Vj g (9) ltered and subsampled to produce the N=4-periodic signals x2 and v2 . The total DWT of x0 is the collection and of signals, v1; : : :; vJ ; xJ . Clearly v1 depends on the highest frequency content of x0, while vJ and xJ depend on BJ = fD y : y = (0; : : :; 0; xJ ); xJ 2 XJ g: (10) the lowest frequency content of x0. Hence, the Decomposition and Reconstruction Algorithms can be implemented 3. Discrete wavelet transforms by analysis and synthesis octave-band lter banks with perTo make the foregoing discussion concrete, let N be a fect reconstruction [26]. power of 2, and set 3.1. Orthonormal wavelet bases Xj = Cj N=2j ; j = 0; : : :; J; We now construct a special orthonormal basis for X0 . Recalling Fact 4 in Section 2.2, it suces to construct orj where 2 N=2J , and Cj N=2 is equipped with the usual thonormal bases for Wj and BJ . j Euclidean inner product. We take Vj = Xj ; j = 1; : : :; J. For each space Cj N=2 , let j denote the N=2j -periodic Kronecker delta, and set ji (n) = j (n ? i). Then the Convention. As in the theory of discrete Fourier transj N= 2 j forms, we do not distinguish between vectors in C and ji, i j= 0; : : :; N=2j ? 1, constitute the standard basis for their doubly-in nite N=2j -periodic extensions. Recall that Cj N=2 = Xj = Vj . We rst construct a basis for Y = 4
4
4
4
3
V1 VJ XJ . Set
wavelet transform is then implemented by the Decomposition Algorithm, and the inverse transform is implemented by the Reconstruction Algorithm. We now show that the complexity of these algorithms is order N. Since gj (n) and hj (n) have only p nonzero terms, the number of multiplications required to compute (11) and (12) for m = 0; : : :; N=2j +1 ? 1 is pN=2j . To implement the Decomposition Algorithm requires that we do this for j = 0; : : :; J ? 1. Hence, the total number of multiplications is always less than 2pN. Further speedups can be achieved by adapting some of the ideas in Rioul and Duhamel [19, Section III]. For example, at each iteration j, the sums in (11) and (12) can be split into sums over n = even and n = odd. This results in the sum of two convolutions. Depending on the relative values of p and log2 (N=2j +1), it may be more ecient to do this with fast Fourier transforms. For small p, Rioul and Duhamel suggest the use of short-length fast running FIR algorithms [16], [25]. The fast wavelet transform can also be expressed as a matrix factorization into sparse matrices analogous to the fast Fourier transform [24]. 3.3. Numerical example To obtain the g(n) mentioned at the beginning of the previous subsection, it is necessary to solve (36){(37) in Section 4. We take p = 4 so that these equations can easily be solved by hand. In the course of the calculations, there are two places where square roots must be taken. Thus, there are four possible solutions, depending on whether positive or negative roots are used. If the negative root is always taken, we obtain
dji = (0; : : :; ji; : : :; 0); j = 1; : : :; J; 4
where ji is in the jth position out of J + 1 positions. Set eJi = (0; : : :; 0; Ji). If we let 4
ji
= D dji and 'Ji = D eJi ; 4
4
(13)
and if we assume that D is unitary, then
f j 0; : : :;
j;N=2j ?1g
and f'J 0; : : :; 'J;N=2J ?1g
are orthonormal bases for Wj and BJ , respectively. The union of all these bases is called a wavelet basis for X0 . The signals ji are called wavelets, and the signals 'Ji are called scaling functions. The connection between ji and 'Ji and the DWT output signals is vj (i) = hx0; jii0 and xJ (i) = hx0 ; 'Jii0 : (14) This is easily seen by writing hx0 ; jii0 = hx0 ; D djii0 = (Dx0 ; dji) = (vj ; ji)j = vj (i): The second equation is obtained similarly. Remarks. (i) The Reconstruction Algorithm shows that 1i = H0 1i , and ji
= G0 Gj ?2Hj?1ji; j = 2; : : :; J:
(15)
Also
'Ji = G0 GJ ?1Ji : (16) (ii) Using (11){(12) and (15){(16), it is an easy calculation to show that for j = 1; : : :; J, j;i+l (n)
=
ji (n ? l2j );
p
and 'J;i+l (n) = 'Ji (n ? l2J ); l = 0; : : :; N=2J ? 1: (18) In particular, it follows that the subspace Wj in (9) (resp. BJ in (10)) is closed under circular shifts of 2j (resp. 2J ) units. 3.2. Fast wavelet transforms Let p be an even integer such that 2 p N=2J . As will be shown in Section 4, one can nd a sequence g(0); : : :; g(p ? 1) such that if n = 0; : : :; p ? 1; (19) gj (n) = g(n); 0; n = p; : : :; N=2j ? 1; and if 4
hj (n) = (?1)n gj (1 ? n); n = 0; : : :; N=2j ? 1; 4
p
g(0) = 1 +p 3 ; g(1) = 3 +p 3 ; 4 2 4 2 (21) p p 1 ? 3 ? 3 3 p ; g(3) = p : g(2) = 4 2 4 2 These are the rst four entries in [9, Table 1, p. 980] and [10, Table 6.1, p. 195], where they are called h(n). If we had always taken the positive square root, we would have obtained g~(n) = g(3 ? n), n = 0; 1; 2; 3. The two remaining solutions are ?g(n) and ?g~(n). Let g be given by (21) and de ne gj and hj by (19) and (20) with N = 64 and J = 4. Note that since gj and hj are real, if x0 is real, so is Dx0 . The signals 1;2 , 2;2, 3;2, 4;2, and '4;1 are shown in Figs. 1{5, respectively. The remaining basis signals can be obtained by circular shifts as indicated in (17) and (18). 4. Construction of Gj and Hj for fast wavelet trans-
l = 0; : : :; N=2j ? 1; (17)
forms
(20)
As we saw in Section 3.2, the key to fast wavelet transwhere j = 0; : : :; J ? 1, then the operators Gj and Hj given forms is nding circular FIR lters gj . Then by (20), hj by (11) and (12), respectively, satisfy (4){(7). The discrete will also be circular FIR. Of course, we also need to show 4
that (4){(7) hold for the corresponding operators Gj and Hj . The plan is to construct Gj from G0 and then to construct Hj from Gj . 4.1. Obtaining Hj from Gj Without loss of generality, we may work with G0 and H0 . To simplify the notation, we drop the subscript 0, and we let M = N=2. Then (11) and (12) become
all u. For each u, let g^o(u) be any complex number such that jg^o(u)j2 = 1 ?jg^e(u)j2. Then (27) holds for all u. If we let h(n) be given by (32), then G and H will satisfy (4){ (7). In the next subsection, we consider the more delicate task of constructing g(n) so that in addition to (27), g is also circular FIR. 4.2. Selection of G0 As before we drop the subscript 0. The purpose of NX ?1 this subsection is to nd a circular FIR solution of (5), or g(n ? 2m)x(n); (22) equivalently, (Gx)(m) = (27). n=0 Using (22), it is easy to see that (5) is equivalent to NX ?1 h(n ? 2m)x(n); (23) (Hx)(m) = NX ?1 n=0 g(n + 2m)g(n) = (m); (33) n =0 where g, h, and x are N-periodic sequences. The purpose of this subsection is to use (22){(23) to express (4){(7) in where is M-periodic. Let (m) denote the left-hand side terms of the lters g and h. Proofs of the following lemma of (33). Clearly, has period M. Also, (M ? m) = (m). and theorem, which are adapted from [9, Section 3.A], are Thus, (m) = 0 for m = 1; : : :; M=2 ? 1 if and only if given in Appendix A. (m) = 0 for m = M=2 + 1; : : :; M ? 1. Therefore, (33) holds if and only if (m) = (m) for m = 0; : : :; M=2. We Lemma 1. De ne the M-periodic sequences ge(k) = now impose the constraint g(n) = 0 for p n N ? 1, g(2k); go (k) = g(2k + 1); he(k) = h(2k), and ho(k) = where p M and p is even. It then follows that (33) can h(2k + 1). Denote their M-length DFTs by g^e(u), g^o(u), be replaced by h^ e (u), and ^ho(u), respectively. If G and H are given by pX ?1 (22){(23), then (4) is equivalent to g(n + 2m)g (n) = (m); m = 0; : : :; M=2: (34) 2 2 n =0 ^ jg^e(u)j + jhe(u)j = 1; (24) ^ ^ g^e (u)^go(u) + he (u)ho(u) = 0; (25) Now, in (34), n+2m (p?1)+2(M=2) = p?1+M N ?1, since we assumed p M. Since there is no wraparound, 2 2 ^ jg^o(u)j + jho(u)j = 1; (26) the sum in (34) will automatically be zero if 2m p. Thus, we may rewrite (34) as and (5){(7) are, respectively, equivalent to 4
4
4
4
?1 (27) pX g(n + 2m)g (n) = (m); m = 0; : : :; p=2 ? 1: (35) (28) n=0 (29) Finally, observe that n + 2m p ? 1 if and only if n g^e (u)^he (u) + g^o (u)^ho(u) = 0: p ? 2m ? 1, and since m p=2 ? 1, p ? 2m ? 1 1. Hence, Furthermore, (24){(26) and (29) together imply (35) becomes j^he(u)j = jg^o(u)j and j^ho(u)j = jg^e(u)j: (30) p?X 2m?1 g(n + 2m)g(n) = (m); m = 0; : : :; p=2 ? 1: Similarly, (27){(29) and (25) together also imply (30). In n=0 addition, (27) and (30) imply (24), (26), and (28). (36) Obviously, (36) does not depend on the value of N M Theorem 2 (Wavelet construction). Let G and H be p. Thus, suppose p is even and that we have a solution to given by (22){(23). If (27) holds, and if h is de ned via (36). Suppose also that N and J are such that 2 p J N=2 if gj is given by (19), then (33) will hold h^ e (u) = g^o (u) and ^ho(u) = ?g^e (u); (31) for gj .. IfClearly, we now de ne hj by (20), Theorem 2 shows that (4){(7) will hold for Gj and Hj if they are de ned by (11) then (4){(7) hold. Furthermore, (31) is equivalent to and (12), respectively. As shown in Section 3.2, we obtain (32) a fast wavelet transform. h(n) = (?1)n g(1 ? n): We now observe that (36) gives us only p=2 nonlinear Clearly, this theorem gives a simple recipe to construct equations in p unknowns. Therefore, we must impose adwavelets. Solutions of (27) are easy to construct. Let g^e (u) ditional consistent constraints on g(0); : : :; g(p ? 1) in orbe arbitrary except for being scaled so that jg^e(u)j2 1 for der to x a solution. The linear constraints suggested by
jg^e(u)j2 + jg^o(u)j2 = 1; j^he(u)j2 + j^ho(u)j2 = 1;
5
Daubechies [9] are
but they are not orthonormal unless D is unitary. Now de ne another basis for X0 by setting
p?1 dk X n dz k n=0 g(n)z z=?1 = 0; k = 0; : : :; p=2 ? 1: (37) Daubechies' original reason for introducing (37) was related to regularizing the behavior of the continuous-time analog of (16) for large J [9, Section 3.B]. However, (37) can also be interpreted as saying that the transfer function P corresponding to g, namely pn?=01 g(n)e ?jn! , as well as its p=2 ? 1 derivatives with respect to !, are zero at ! = . In other words, this is a low-pass transfer function with a maximally at stopband. Clearly, the N-length DFT of g(n), denoted by g^(u), is obtained by sampling this transfer function at ! = 2u=N. If h(n) is then de ned by (32), it is easy to see that the corresponding transfer function is high-pass. Theorem 3 (Daubechies). The constraints (37) are consistent with (36). Real-valued solutions g(n) exist, some of which are given in [10, Table 6.1], where they are called h(n). The proof of this result, which is e ectively contained in [9, Section 4.B], is sketched in Appendix A. The theorem statement only guarantees existence of a solution to (36){ (37). In fact, the proof is constructive, obtaining solutions by spectral factorization of certain polynomials; however, this is not the only way to obtain a solution in a given instance. One can attempt a direct solution by hand if the system of equations is small (e.g., p = 4 as in Section 3.3) or by numerical techniques such as Newton's method otherwise. Note that if g(n) is real, then so is h(n). Hence, D can also be viewed as a mapping from IRN to IRN .
= D?1 dji and 'eJi = D?1 eJi : 4
ji
e
4
These signals and the signals in (13) form a biorthogonal set in the sense that h ji ; ej i i0 = hD dji; D?1dj i i0 = (dji; DD?1 dj i ) = ii jj , h ji; 'eJi i0 = h'Ji ; eji i0 = 0, and h'Ji ; 'eJi i0 = ii . From (8) and the de nition of dji and eJi in Section 3.1, we can write 0 0
0 0
0
0 0
0
0
J
N=X 2j ?1
j =1
i=0
X
Dx0 =
0
0
0
vj (i)dji +
N=X 2J ?1 i=0
xJ (i)eJi :
Recalling (14), and noting that x0 = D?1 (Dx0 ), we have the following generalization of (1): x0 =
J
N=X 2j ?1
j =1
i=0
X
hx0 ; jii0
ji
e
+
N=X 2J ?1 i=0
hx0; 'Ji i0 'eJi :
e Let Dx 0 denote the output of the Decomposition Algorithm when Gj is replaced by Aj and Hj is replaced by Cj. e Then De = D?1 and De ?1 = D . Since x0 = De ?1 (Dx 0),
x0 =
J
N=X 2j ?1
j =1
i=0
X
hx0 ; ejii0
ji
+
N=X 2J ?1 i=0
hx0; 'eJi i0 'Ji :
Note also that ji = D dji = D (D eji ) = (D D) eji , and so eji = (D D)?1 ji. This further implies eji = e ?1 e . Of course, ' e e D) (De D) ji Ji and 'eJi ji and ji = (D are similarly related. We say that f ji; 'Ji g and f eji; 'eJig 5. Biorthogonal wavelets are dual wavelet bases (cf. [10, Eq. (3.3.16)]). We now summarize the generalizations of Section 4.1. The wavelets shown in Figs. 1{5 are obviously asymFirst we must augment (22) and (23) with (recall M = metric. Since the Decomposition and Reconstruction AlN=2) gorithms correspond to exact reconstruction lter banks in which the analysis and synthesis lters use the same coeMX ?1 cients, this asymmetry is to be expected [22]. The diculty a(n ? 2m)x(m); (Ax)(n) = arises because we were too restrictive in Section 2.1 when m =0 we required T T = I and TT = I. We can relax these MX ?1 conditions by requiring only that T be invertible, where c(n ? 2m)v(m); (Cv)(n) = T ?1 has the form T ?1 = [A C ]. Then the generalizam=0 tions of (4){(7) are AG + CH = I, GA = I, HC = I, and GC = 0, respectively. We also need to add HA = 0. where the complex conjugate is included only for conveThe Decomposition Algorithm, which uses only Gj and nience. After reading the proof of Lemma 1 in Appendix A, Hj , stays the same, and its output is still called Dx0. The it is easy to see that in the present context, AG + CH = I Reconstruction Algorithm, which uses only Gj and Hj , is equivalent to still computes D , but now D 6= D?1 . However, we can (38) g^e(u)^ae (u) + ^he (u)^ce (u) = 1; compute D?1 if we modify the Reconstruction Algorithm by replacing Gj by Aj and Hj by Cj , assuming that Aj Gj + ^ (39) g^e (u)^ao(u) + he (u)^co (u) = 0; Cj Hj = I for each j. ^ g^o(u)^ao (u) + ho (u)^co (u) = 1; (40) If D is invertible, then so is D . In this case, the signals N ^ j g^o (u)^ae(u) + ho (u)^ce (u) = 0; (41) ji and 'Ji de ned in (13) still form a basis for X0 = C , 6
where (p; p~; m) = minfp ? 2m ? 1; p~ ? 1g. (If p~ = 1 or 2, we have only the rst equation.) Clearly, these conditions will hold for all N M maxfp; p~; 2g. 5.1. Symmetric, biorthogonal, FIR wavelets Let ( (maxfp; p~g ? 1)=2; p; p~ both odd; K = maxfp; p~g=2 ? 1; p; p~ both even: We view both g(n) and a(n) as being of length 2K + 1 in the rst case, and of length 2K + 2 in the second case. In other words, we pad the tail of the shorter sequence with zeros. We now require that g and a be symmetric. For p and p~ both odd, this means that g(K + n) = g(K ? n) and a(K + n) = a(K ? n) for n = 1; : : :; K. For p and p~ both even, this means that g(K +n) = g([K +1] ? n) and a(K + n) = a([K+1] ? n) for n = 1; : : :; K+1. Note that since we zero-padded the tail of the shorter sequence, by imposing symmetry, we are forcing zeros in the head as well. To pin down choices further for g and a, we impose conditions analogous to (37). With slight abuse of notation, let
and GA = I, HC = I, GC = 0, HA = 0 are, respectively, equivalent to g^e (u)^ae (u) + g^o(u)^ao (u) = 1; (42) ^he (u)^ce (u) + h^ o (u)^co (u) = 1; (43) g^e (u)^ce (u) + g^o (u)^co (u) = 0; (44) ^he (u)^ae (u) + h^ o (u)^ao (u) = 0: (45) Suppose that g and a are such that (42) holds. Taking our hint from (31), let h and c be de ned via h^ e (u) = a^o (u) and h^ o(u) = ?a^e (u); (46)
4
4
and
c^e (u) = g^o(u) and c^o (u) = ?g^e (u): (47) Then (38){(41) and (43){(45) also hold. Note that (46) and (47) are equivalent to h(n) = (?1)n a(1 ? n) and c(n) = (?1)n g(1 ? n); (48) respectively. Thus, we now have a simple recipe to construct biorthogonal wavelets. If g^e (u) and g^o (u) are arbitrary except for the constraint that g^e (u)^go(u) 6= 0 for all u, then we can always nd choices of a^e (u) and a^o(u) that satisfy (42). If h and c are then de ned by (48), and if we can do this for j = 0; : : :; J ? 1, then we have a system of biorthogonal wavelets. In order that the algorithms for computing D, D , and ? 1 D be fast, we would like to construct compactly supported biorthogonal wavelets based on two sequences, to be chosen below, say g(0); : : :; g(p ? 1) and a(0); : : :; a(~p ? 1), where we assume that 1 minfp; p~g and maxfp; p~; 2g N=2J . We do not require that p or p~ be even. We de ne the N=2j -periodic sequence gj by (19), and de ne aj similarly. We use (48) to de ne hj in terms of aj and cj in terms of gj . To be sure this will all work correctly, we must generalize the derivation leading from (33) to (36). Now, GA = I and (42) are both equivalent to NX ?1
G(z) = 4
pX ?1 n=0
g(n)z n and A(z) = 4
pX ~?1 n=0
a(n)z n : (51)
We require that the derivatives with respect to z satisfy G(k)(?1) = 0; k = 0; : : :; ` ? 1; (52) A(k) (?1) = 0; k = 0; : : :; `~ ? 1; where 1 ` p and 1 `~ p~. Note that these conditions do not uniquely specify a solution of (49); if g(n) and a(n) satisfy (50) and (52), then so do g(n)= and a(n) for any nonzero constant . We therefore impose the normalization condition p~?1 X p (53) a(n) = 2:
p
n=0
The choice of 2 is only for convenience; any nonzero con(49) stant would do. n=0 Example. Let p = 6 and p~ = 4. Then K = 2. By imposwhere g and a are regarded as N-periodic and is M = ing symmetry, only a(3) = a(2), g(5) = g(0), g(4) = g(1), N=2-periodic. Let (m) denote the left-hand side of (49). and g(3) = g(2) are nonzero. Substituting these into (50) Observe that (m) = (m) for m = 0; : : :; M ? 1 if and only yields the two independent equations, g(2)a(2) = 12 and if (m) = (m) for m = 0; : : :; M=2 and (M ? m) = 0 g(0) = ?g(1). Take ` = 2 and put no constraint on A(z). for m = 1; : : :; M=2 ? 1. If we now impose the conditions Since g is symmetric, G(0)(?1) = 0 already. G0 (?1) = 0 g(n) = 0 for n = p; : : :; N ? 1 and a(n) = 0 for n = implies p g(2) = 8g(1). The normalization (53) implies p~; : : :; N ? 1, and proceed as in Section 4.2, we see that a(2) = p2=2. Thus, g(2) = p2=2, g(1) = p2=16, and (49) is equivalent to the pair of equations, g(0) = ? 2=16. The remaining coecients are pfound from symmetry. If we divide all our coecients by 2, we ob(p; p~;m) X tain the coecients in the second row of Table 8.2 in [10, g(n + 2m)a(n) = (m); 0 m < p=2; p. 277]. For N = 64 and J = 3, plots of e3;2, 'e3;1, 3;2, n=0 (50) and '3;1 are shown in Figs. 6{9, respectively. (~pX ;p;m) When we imposed (52), there was no guarantee that g(n)a(n + 2m) = 0; 1 m < p~=2; these constraints were consistent with (42). That this is n=0 g(n + 2m)a(n) = (m);
7
indeed the case essentially follows from [10, Section 8.3.4]. Recall that (36) was derived from (33), and (33) is We brie y sketch the argument in Appendix A. equivalent to (27). Letting g^(u) denote the N-length DFT of g(n), and using the analog of (55), we nd that (27) is Appendix A equivalent to Proof of Lemma 1. To show that (4) is equivalent to jg^(u)j2 + jg^(u + N=2)j2 = 2: (57) (24){(26), rst note that (22){(23) imply that (4) is equivalent to Now observe that g^(u) = G(e j! ) when ! = ?2u=N and u is an integer. If (56) holds, then MX ?1 g(n ? 2m)g( ? 2m)+ h (n ? 2m)h( ? 2m) = (n ? ); jG(e j! )j2 = [cos2(!=2)]p=2 jF(e j! )j2: m=0 (54) Next, if F(z) has real coecients, then jF(e j! )j2 also has where denotes the N-periodic Kronecker delta. By con- real coecients and is real-valued. Hence, we can write sidering (54) when n and are both even, it is easy to jF(e j! )j2 = F(e j! )F(e ?j! ) as a polynomial of degree P ?1 see that (54) implies mM=0 ge(k ? m)ge (?m) + he(k ? p=2 ? 1 in cos !, or equivalently, in sin2 (!=2); i.e., we may m)he (?m) = (k), where here is M-periodic. Tak- write jF(e j! )j2 = P(sin2 (!=2)), where P is a polynomial ing an M-length DFT and changing ?u to u yields (24). of degree p=2 ? 1. So, if (56) holds for some F with real Similarly, taking n odd and even in (54) leads to (25). coecients, and if (57) holds (which is what we are trying Taking n and both odd leads to (26). To show that to prove), then (5){(7) are, respectively, equivalent to (27){(29), proceed y p=2P(1 ? y) + (1 ? y)p=2 P(y) = 2; (58) as follows. We treat only (5) as the other two are similar. As noted in Section 4.2, (5) is equivalent to (33). where y = cos2 (!=2), ! = ?2u=N, and u is an integer. Breaking the sum in (33) into even and odd n yields PM ?1 Conversely, the polynomial k=0 ge(k+m)ge (k)+go (k+m)go (k) = (m). Taking an M-length DFT of this equation gives (27). The remaining p=X 2?1 p=2 ? 1 + k statements of the theorem are easily veri ed. 2 yk (59) P(y) = 2 k k =0 Proof of Theorem 2. First note that (31) trivially implies (25) and (29). Next, (31) implies (30), and by hy- satis es (58).1 Furthermore, P(y) is nonnegative for 0 pothesis, (27) holds. Hence, (24), (26), and (28) hold by y 1, and by Riesz's lemma, there exists a (nonunique) Lemma 1. Since (24){(29) now hold, Lemma 1 further im- polynomial Q of degree p=2 ? 1 with real coecients that plies that (4){(7) hold. To show that (31) and (32) are satis es jQ(e j! )j2 = P(sin2(!=2)). Thus, if we take the Q equivalent, substitute (31) into the easily veri ed formula, of Daubechies2 and obtain the g(n) from ^h(u) = ^he (u) + h^ o (u)e ?j2u=N ; G(z) = [ 21 (1 + z)]p=2 Q(z); (55) where h^ is the N-length DFT of h. We obtain h^ (u) = and set ! = ?2u=N and y = cos2 (!=2), then g^o (u) ? g^e (u)e ?j2u=N , from which (32) follows. The fact jg^(u)j2 + jg^(u + N=2)j2 that (32) implies (31) follows by separately considering (32) = jG(e j! )j2 + jG(e j(!?) )j2 with n replaced by 2n + 1 and by 2n. 2 = y p=2P(1 ? y) + (1 ? y)p=2 P(y) Proof of Theorem 3. The proof is self-contained ex= 2: 2 cept for the reference to Riesz's lemma (see e.g., [9, p. 972, Lemma 4.2] or [10, Lemma 6.1.3]), which is a result about Proof of Consistency of (52). We treat the case in the spectral factorization of polynomials that can be read which p and p~ are both even, the other case being similar. independently of the rest of [9] or [10]. First note that just as (57) was established above, (42) is Let G(z) be as in (51). Since (37) is equivalent to saying equivalent to that (1 + z)p=2 is a factor of G(z), g^(u)^a(u) + g^(u + N=2)^a(u + N=2) = 2: (56) G(z) = [ 21 (1 + z)]p=2F(z) 1 ?p=2 p=2 4
4
4
If P (y) satis es (58), P (y) = (1 ? y) [2 ? y P (1 ? y)]. Substituting the Taylor series of (1 ? y)?p=2 yields (59) plus higherorder terms. Since P (y) is a polynomial of degree less than p=2, the higher-order terms must sum to zero. Eq. (59) is the unique lowestdegree solution of (58). Higher-degree solutions also exist [10, p. 171].
for some polynomial F of degree p=2 ? 1. The factor of 1=2 is included for convenience in the derivation below. Observe that the coecients of G are all real if and only if the coecients of F are all real. We now require that 2 In general there are p=2 choices for Q(z ) (and hence G(z )). The g(n) be real. Thus, it suces to prove the existence of a coecients (n) (called h(n) in [9], [10]) given in [9, p. 980, Table 1], polynomial F of degree p=2 ? 1 with real coecients such [10, p. 195, gTable 6.1] correspond to taking Q(z) to have all its roots that if G(z) is de ned by (56), then (36) holds. inside the unit circle; i.e., Q(z) is a minimum-phase lter. 8
If we set ! = ?2u=N, this can be written as G(e j! )A(e j! ) + G(e j(!?) )A(e j(!?) ) = 2:
[5] G. Caire, R. L. Grossman and H. V. Poor, \Wavelet transforms associated with nite cyclic groups", IEEE Trans. Inform. Theory, Vol. 39, No. 4, July 1993, pp. 1157{1166. [6] C. K. Chui, An Introduction to Wavelets, Academic Press, New York, 1992. [7] A. Cohen, I. Daubechies and J.-C. Feauveau, \Biorthogonal bases of compactly supported wavelets", Commun. Pure Appl. Math., Vol. 45, 1992, pp. 485{560. [8] R. R. Coifman and M. V. Wickerhauser, \Best-adapted wavelet packet bases", preprint, February 1990. [9] I. Daubechies, \Orthonormal bases of compactly supported wavelets", Commun. Pure Appl. Math., Vol. 41, November 1988, pp. 909{996. [10] I. Daubechies, Ten Lectures on Wavelets, SIAM, Philadelphia, PA, 1992. [11] I. Gohberg and S. Goldberg, Basic Operator Theory, Birkhauser, Boston, 1980. [12] C. E. Heil and D. F. Walnut, \Continuous and discrete wavelet transforms", SIAM Rev., Vol. 31, No. 4, December 1989, pp. 628{666. [13] F. Hlawatsch and G. F. Boudreaux-Bartels, \Linear and quadratic time-frequency signal representations", IEEE Signal Process. Mag., Vol. 9, No. 27, April 1992, pp. 21{67. [14] K. Ho man and R. Kunze, Linear Algebra, 2nd Edition, Prentice-Hall, Englewood Cli s, NJ, 1971. [15] S. Mallat, \Multiresolution approximations and wavelet orthonormal bases of L2 (IR)", Trans. Amer. Math. Soc., Vol. 135, No. 1, September 1989, pp. 69{88. [16] Z. J. Mou and P. Duhamel, \Short-length FIR lters and their use in fast nonrecursive ltering", IEEE Trans. Signal Process., Vol. 39, June 1991, pp. 1322{1332. [17] W. H. Press, S. A. Teukolsky, W. T. Vetterling and B. P. Flannery, Numerical Recipes in C: The Art of Scienti c Computing, 2nd Edition, Cambridge Univ. Press, New York, 1992. [18] O. Rioul, \A discrete-time multiresolution theory", IEEE Trans. Signal Process., Vol. 41, No. 8, August 1993, pp. 2591{ 2606. [19] O. Rioul and P. Duhamel, \Fast wavelet algorithms for discrete and continuouswavelet transforms", IEEE Trans. Inform. Theory, Vol. 38, No. 2, March 1992, pp. 569{586. [20] O. Rioul and M. Vetterli, \Wavelets and signal processing", IEEE Signal Process. Mag., Vol. 8, No. 4, October 1991, pp. 14{ 38. [21] M. J. Shensa, \The discrete wavelet transform: Wedding the a trous and Mallat algorithms", IEEE Trans. Signal Process., Vol. 40, No. 10, October 1992, pp. 2464{2482. [22] M. J. T. Smith and T. P. Barnwell III, \Exact reconstruction techniques for tree-structured subband coders", IEEE Trans. Acoust. Speech, Signal Process., Vol. ASSP-34, No. 3, June 1986, pp. 434{441. [23] G. Strang, \Wavelets and dilation equations: A brief introduction", SIAM Rev., Vol. 31, No. 4, December 1989, pp. 614{627. [24] G. Strang, \Wavelet transforms versus Fourier transforms", Bull. Amer. Math. Soc., Vol. 28, No. 2, April 1993, pp. 288{ 305. [25] M. Vetterli, \Running FIR and IIR ltering using multirate lter banks", IEEE Trans. Acoust. Speech, Signal Process., Vol. ASSP-36, May 1988, pp. 730{738. [26] M. Vetterli and C. Herley, \Wavelets and lter banks: Theory and design", IEEE Trans. Signal Process., Vol. 40, No. 9, September 1992, pp. 2207{2232.
(60)
Now, the symmetry of g(n) and a(n) implies G(e j! ) = e j!(K +1=2) cos(!=2) q(cos !); A(e j! ) = e j!(K +1=2) cos(!=2) r(cos !);
(61)
where q and r are polynomials of degree less than or equal to K. If (52) also holds, then G(z) = (1 + z)` G0(z) and A(z) = (1 + z)`~A0(z) (62) for some polynomials G0 (z) and A0 (z). It can be shown that (61){(62) imply that G(e j! ) = e j!(K +1=2) cos+1 (!=2) q0(cos !); A(e j! ) = e j!(K +1=2) cos~+1 (!=2) r0 (cos !); where q0 and r0 are polynomials, and where = ` if ` is even and = ` ? 1 if ` is odd and similarly for ~. Substituting this into (60), we obtain cos2L(!=2) q0 (cos !) r0 (cos !) + sin2L (!=2) q0(? cos !) r0(? cos !) = 2; where L = ( + ~ + 2)=2 is a positive integer. Since cos ! = 1 ? 2 sin2(!=2), we can write q0(cos !)r0(cos !) = P(sin2 (!=2)) for some polynomial P. Letting y = sin2 (!=2), we have (1 ? y)L P(y) + yL P(1 ? y) = 2: As shown above in the proof of Theorem 3, there are polynomials P that satisfy this equation. By factoring P, feasible q0 and r0 can be found. Note that the choice of factors imposes a normalization condition, though not necessarily the one in (53). 2
Acknowledgements
This work was supported in part by the Air Force Oce of Scienti c Research under Grants AFOSR{90{0181 and F49620{92{J{0305.
References
[1] B. Alpert, G. Beylkin, R. Coifman and V. Rokhlin, \Waveletlike bases for the fast solution of second-kind integral equations", SIAM J. Sci. Comput., Vol. 14, No. 1, January 1993, pp. 159{184. [2] G. Beylkin, \On the representation of operators in bases of compactly supported wavelets", preprint. [3] G. Beylkin, R. Coifman and V. Rokhlin, \Fast wavelet transforms and numerical algorithms I", Comm. Pure Appl. Math., Vol. 44, 1991, pp. 141{183. [4] P. J. Burt and E. H. Adelson, \The Laplacian Pyramid as a compact image code", IEEE Trans. Commun., Vol. COM-31, April 1983, pp. 532{540.
9
1 0.9 0.8 0.7 0.6 0.5 0.4
0.5
0.3 0.4
0.2 0.1
0.3
0 -0.1
0.2
-0.2 0.1
-0.3 -0.4
0
-0.5 0
10
20
30
40
50
-0.1
60
-0.2
Fig. 1. Wavelet 1;2(n) of Section 3.3.
-0.3 -0.4 0
10
20
30
40
50
60
Fig. 4. Wavelet 4;2(n) of Section 3.3.
0.8 0.7 0.6 0.5 0.4 0.3 0.2 0.1 0 -0.1 -0.2 -0.3 -0.4 -0.5 -0.6 0
10
20
30
40
50
60
Fig. 2. Wavelet 2;2(n) of Section 3.3.
0.40 0.35 0.30 0.25
0.6
0.20
0.5
0.15
0.4
0.10
0.3 0.05 0.2 0.1
0
0
-0.05
-0.1
-0.10
-0.2
0
10
20
30
40
50
60
-0.3
Fig. 5. Scaling function '4;1(n) of Section 3.3.
-0.4 -0.5 -0.6 0
10
20
30
40
50
60
Fig. 3. Wavelet 3;2(n) of Section 3.3. 10
0.4 0.3 0.2 0.1 0 -0.1 -0.2 -0.3 -0.4 0
10
20
30
40
50
60
Fig. 6. Biorthogonal wavelet e3;2(n) of Section 5.
0.45
0.40
0.40
0.35
0.35 0.30
0.30
0.25
0.25 0.20
0.20 0.15 0.15
0.10
0.10
0.05 0
0.05
-0.05
0
-0.10 0
10
20
30
40
50
0
60
10
20
30
40
50
60
Fig. 7. Biorthogonal scaling function 'e3;1(n) of Section 5. Fig. 9. Biorthogonal scaling function '3;1(n) of Section 5.
0.5 0.4 0.3 0.2 0.1 0 -0.1 -0.2 -0.3 -0.4 -0.5 0
10
20
30
40
50
60
Fig. 8. Biorthogonal wavelet 3;2(n) of Section 5. 11
Wavelet transforms for discrete-time periodic signals
John A. Gubner and Wei-Bin Chang
Department of Electrical and Computer Engineering, University of Wisconsin, Madison, WI 53706, USA
Abstract. Wavelet transforms for discrete-time periodic signals are developed. In this nite-dimensional context, key ideas from the continuous-time papers of Daubechies and of Cohen, Daubechies, and Feauveau are isolated to give a concise, rigorous derivation of the discrete-time periodic analogs of orthonormal and symmetric biorthogonal bases of compactly supported wavelets. These discrete-time periodic wavelets are expressed in terms of circular FIR lters, and thus lead to fast wavelet transforms whose complexity is order N . Keywords: Orthonormal wavelets; Biorthogonal wavelets; Symmetric biorthogonal wavelets; Pyramid algorithms; Fast wavelet transforms
1. Introduction
tion 2 we review the properties of abstract pyramid-type algorithms [4], which are the key to fast wavelet transforms. Section 3 specializes these abstract pyramid algorithms to discrete-time periodic signals, de nes the discrete wavelet transform, and discusses fast wavelet transforms. Section 4 derives the discrete-time periodic analogs of Daubechies' orthonormal bases of compactly supported wavelets [9]. Section 5 derives the discrete-time periodic analogs of the symmetric, biorthogonal bases of compactly supported wavelets in [7]. 1.1. Sinusoids and wavelets Let X0 denote the space of discrete-time signals x0 (n) with period N, where N is a power of 2. Before de ning wavelets, recall that from the theory of discrete Fourier transforms, every signal x0 2 X0 can be expressed in terms of N complex sinusoids as
Although there is a large and growing literature on wavelets, there is, to our knowledge, no reference that focuses on the development of wavelets for discrete-time periodic signals, and that includes a derivation of Daubechies' FIR lter coecients in this nite-dimensional context. This paper is an attempt to ll this gap, which is making its presence felt as researchers use discrete-time periodic wavelets in practical applications ranging from image compression [17, Section 13.10] to numerical analysis [1]{[3], [8]. Because discrete-time periodic signals constitute a nitedimensional vector space, our presentation is based primarily on concepts from linear algebra and discrete Fourier transforms, though some details in Appendix A require results about the factorization of polynomials. Hence, we are able to give a concise, rigorous exposition of the key ideas in [9] and [7] that pertain to the discrete-time periodic analogs of compactly supported orthonormal and symmetric biorthogonal wavelets. Our presentation is aimed at the non-wavelet-specialist who is not familiar with the more technically demanding continuous-time papers [9] and [7]. The relationship between the continuous wavelet transform and the discrete wavelet transform for discrete-time aperiodic signals has been studied by Shensa [21]. Multiresolution analysis for continuous-time signals was developed by Mallat [15], and multiresolution analysis for discrete-time aperiodic signals was developed by Rioul [18]. Wavelet transforms for nite-length sequences have been developed by Caire et al. [5], and some of their equations can be found in our paper. However, [5] focuses on signals taking values in nite elds and does not consider the consistency of Daubechies' constraint (37). Also, [5] does not consider biorthogonal wavelets for discrete-time periodic signals as we do in Section 5. The connection between wavelets and lter banks has been well documented by Vetterli and Herley [26]. Recent expository treatments of wavelets include [6], [10], [12], [13], [20], [23], and [24]. The paper is organized as follows. In the remainder of this section, we introduce discrete-time periodic wavelets and compare them to discrete-time sinusoids. In Sec-
x0(n) = where
j2kn=N
NX ?1 k=0
hx0; Ek iEk (n); NX ?1
Ek (n) = e p ; hx0 ; Eki = x0 (n)Ek (n); N n=0 and the overbar denotes the complex conjugate. The rea?1 son the decomposition works is that the signals fEk gNk=0 form an orthonormal basis for X0 . While direct computation of the Fourier coecients hx0 ; Eki for k = 0; : : :; N ? 1 would require N 2 multiplications, the fast Fourier transform can compute these coecients with order N log2 N multiplications. In this paper we consider orthonormal \wavelet" bases for X0 . To create a wavelet basis of dimension N, we rst choose the number of scales or resolution levels, denoted by J, where N=2J 2. The basis is then partitioned according to resolution level as follows. At level j, j = 1; : : :; J, there are N=2j wavelet signals denoted by j;0; : : :; j;N=2j ?1 . As this accounts for only N ? N=2J signals, the basis is completed by adding N=2J scaling functions, 'J;0 ; : : :; 'J;N=2J ?1 . Since these signals are chosen 4
1
4
to form an orthonormal basis, every x0 2 X0 has the ex- uct on U V by taking pansion u ; u~ = hu; u~i + (v; v~): v v~ 2j ?1 J N=X X hx0 ; jii ji (n) x0(n) = Now, if X is a third inner-product space, and if G: X ! U j =1 i=0 (1) and H: X ! V are linear transformations, we de ne the N=X 2J ?1 mapping T: X ! U V by hx0; 'Ji i'Ji (n): + G x = Gx : Tx = H (3) i=0 Hx The numbers hx0; jii and hx0 ; 'Jii are called wavelet co- Next, if the adjoints G : U ! X and H : V ! X exist, ecients, and the fast wavelet transform is an algorithm then the adjoint T : U V ! X also exists and is given by for computing them with order N multiplications. Ex u u amples of wavelets for the case N = 64 and J = 4 are T v = [G H ] v = Gu + H v: shown in Figs. 1{5. It will be shown later that for each j, ji(n) = j 0 (n ? i2j ) for i = 0; : : :; N=2j ? 1. In Our interest is in isometries; i.e., operators that preserve other words, the ji are just circular translates of j 0 inner products, or equivalently, operators for which T T = by multiples of 2j . Similarly, 'Ji (n) = 'J 0(n ? i2J ) for I. For such operators, x can obviously be recovered from i = 0; : : :; N=2J ? 1. We also point out that a quick glance Tx. at the gures might lead one to wonder if in analogy with Remark. If T T = I and T is onto, then TT = I as continuous time, well [11, p. 185]. Such an operator is said to be unitary. Of ( p course, if dim X = dim U V < 1, then T T = I always 2 j +1;2(2n); j +1;2(2n) 6= 0; ? (2) implies TT = I, since in this case T is 1-to-1 if and only j;2(n) = 0; otherwise: if T is onto [14, p. 81]. However, a close inspection of Figs. 1 and 2 reveals that For future reference, note that since T has the form in this cannot be true since 2;2 has 10 nonzero components, (3), T T = I is equivalent to while 1;2 has only 4 nonzero components. To understand G G + H H = I; (4) more fully what is happening, it is necessary to consider the limiting continuous-time case [9], [10]. and TT = I is equivalent to the following three equations: We now compare and contrast the sinusoidal basis funcGG = I; (5) tions Ek with thepwavelets ji and 'Ji . The constant func HH = I; (6) tion, E0(n) 1= N, exhibits no variation at all, while the signals 'Ji (n) exhibit the least variation of all the wavelet GH = 0: (7) signals. At the other extreme, the signal EN ?1(n) exhibits To summarize, if we \decompose" x with the most variation of all the sinusoids Ek (n), while the signals 1;i(n) exhibit the most variation of all the wavelet u = Gx; v = Hx; signals. One important di erence between the wavelets in and, if (4) holds, then we can reconstruct x from u and v Figs. 1{5 and the sinusoids Ek is that the ji(n) are zero with for most values of n. For example, in Fig. 1, we see that x = Gu + H v: 1;2 (n) is nonzero only for n = 2; 3; 4; 5. Thus, looking at wavelet coecients can help localize high-variation features 2.2. The general algorithms in time rather precisely. Slower-varying signals cannot be Let X0 ; X1; : : :; XJ and V1 ; : : :; VJ be two sequences localized as well. of inner-product spaces. Denote the inner product on Xj by h; ij and the inner product on Vj by (; )j . For 2. Pyramid algorithms j = 0; : : :; J ? 1, let Gj : Xj ! Xj +1 and let Hj : Xj ! Vj +1 The basic building block for pyramid-type algorithms be linear operators (for example, see the beginning of Secis a linear transformation from an input space X into a tion 3). Consider the following procedure. product space, say U V , as discussed in Section 2.1. In Section 2.2 we iterate a sequence of these transformations Decomposition Algorithm to obtain a Decomposition Algorithm and a Reconstruction Let x0 2 X0 be given. Algorithm. for j = 0 to J ? 1 xj +1 = Gj xj 2.1. Linear algebra preliminaries vj +1 = Hj xj Let U and V be inner-product spaces with respective next j inner products h; i and (; ). This induces an inner prod- end 4
2
If (4) holds for each pair Gj and Hj , then x0 can be the N-length discrete Fourier transform (DFT) of a se?1 x(n)e ?j2nu=N , recovered from v1 ; : : :; vJ ; xJ by the following procedure. quence x 2 Cj N is de ned by x^(u) = PNn=0 for all integers u. Reconstruction Algorithm De nition (Discrete wavelet transform). For j = Let v1 2 V1 ; : : :; vJ 2 VJ ; xJ 2 XJ be given. 0; : : :; J ? 1, if the linear operators Gj : Xj ! Xj +1 and for j = J ? 1 to 0 H : Xj ! Vj +1 have the form j xj = Gj xj +1 + Hj vj +1 next j N=X 2j ?1 end gj (n ? 2m)x(n); (11) (Gj x)(m) = n=0 If we let N=X 2j ?1 hj (n ? 2m)x(n); (12) (Hj x)(m) = Y = V1 V J XJ ; 4
n=0
and if v1; : : :; vJ , and xJ are generated by the Decomposi- where g , h , and x are N=2j -periodic sequences, and where j j tion Algorithm, then setting gj and hj are such that (4){(7) hold for Gj and Hj , then Dx0 = (v1; : : :; vJ ; xJ ) (8) Dx0 de ned by using these operators in the Decomposition Algorithm is called a discrete wavelet transform (DWT) of de nes a mapping D: X0 ! Y . We equip Y with the in- x0. duced inner product, denoted by (; ). The induced norm We rst note that (11) and (12) are just subsampled is denoted by k k. Note that we are now using row-vector circular convolutions of x(n) with gj (?n) and hj (?n), renotation for product-space vectors. spectively. Next, it can be seen from Section 4.1 and from With the preceding setup, the following facts are easily (57) in Appendix A, that if Gj and Hj satisfy (4){(7), then veri ed. gj and hj are discrete-time periodic analogs of the conju(1) Assuming only that Gj and Hj exist, the vector x~0 gate quadrature lters in [22]. In Section 4 we construct generated by applying the Reconstruction Algorithm gj and hj under two additional constraints. The rst is to an arbitrary vector y~ = (~v1 ; : : :; v~J ; x~J ) 2 Y , is that gj (n) and hj (n) be circular FIR lters, i.e., they are equal to D y~; i.e., the Reconstruction Algorithm im- zero for most values of n in one period. The second is that plements D . Daubechies' constraint given below in (37) hold, forcing gj (2) If (4) holds for each pair Gj and Hj , then D D = I; to be a low-pass lter with maximally at stopband. It i.e., D is an isometryPand preserves norms, and by (8), then follows that on account of (7), hj will be a high-pass kx0 k20 = kDx0 k2 = Jj=1 kvj k2j + kxJ k2J . lter. When gj is a low-pass lter and hj is a high-pass (3) If (5){(7) hold for each pair Gj and Hj , then DD = I. lter, we see that the Decomposition Algorithm rst takes If (4) also holds, then D is unitary. the N-periodic input x0 and computes x1 and v1 , which (4) If D is unitary, then we can write X0 as a direct sum are both N=2-periodic. Clearly, x1 is the result of low-pass of orthogonal subspaces, X0 = W1 WJ BJ , ltering x0 and then subsampling. Similarly, v1 is the rewhere sult of high-pass ltering x0 and then subsampling. The next time through the loop, x1 is low-pass and high-pass Wj = fD y : y = (0; : : :; 0; vj ; 0; : : :; 0); vj 2 Vj g (9) ltered and subsampled to produce the N=4-periodic signals x2 and v2 . The total DWT of x0 is the collection and of signals, v1; : : :; vJ ; xJ . Clearly v1 depends on the highest frequency content of x0, while vJ and xJ depend on BJ = fD y : y = (0; : : :; 0; xJ ); xJ 2 XJ g: (10) the lowest frequency content of x0. Hence, the Decomposition and Reconstruction Algorithms can be implemented 3. Discrete wavelet transforms by analysis and synthesis octave-band lter banks with perTo make the foregoing discussion concrete, let N be a fect reconstruction [26]. power of 2, and set 3.1. Orthonormal wavelet bases Xj = Cj N=2j ; j = 0; : : :; J; We now construct a special orthonormal basis for X0 . Recalling Fact 4 in Section 2.2, it suces to construct orj where 2 N=2J , and Cj N=2 is equipped with the usual thonormal bases for Wj and BJ . j Euclidean inner product. We take Vj = Xj ; j = 1; : : :; J. For each space Cj N=2 , let j denote the N=2j -periodic Kronecker delta, and set ji (n) = j (n ? i). Then the Convention. As in the theory of discrete Fourier transj N= 2 j forms, we do not distinguish between vectors in C and ji, i j= 0; : : :; N=2j ? 1, constitute the standard basis for their doubly-in nite N=2j -periodic extensions. Recall that Cj N=2 = Xj = Vj . We rst construct a basis for Y = 4
4
4
4
3
V1 VJ XJ . Set
wavelet transform is then implemented by the Decomposition Algorithm, and the inverse transform is implemented by the Reconstruction Algorithm. We now show that the complexity of these algorithms is order N. Since gj (n) and hj (n) have only p nonzero terms, the number of multiplications required to compute (11) and (12) for m = 0; : : :; N=2j +1 ? 1 is pN=2j . To implement the Decomposition Algorithm requires that we do this for j = 0; : : :; J ? 1. Hence, the total number of multiplications is always less than 2pN. Further speedups can be achieved by adapting some of the ideas in Rioul and Duhamel [19, Section III]. For example, at each iteration j, the sums in (11) and (12) can be split into sums over n = even and n = odd. This results in the sum of two convolutions. Depending on the relative values of p and log2 (N=2j +1), it may be more ecient to do this with fast Fourier transforms. For small p, Rioul and Duhamel suggest the use of short-length fast running FIR algorithms [16], [25]. The fast wavelet transform can also be expressed as a matrix factorization into sparse matrices analogous to the fast Fourier transform [24]. 3.3. Numerical example To obtain the g(n) mentioned at the beginning of the previous subsection, it is necessary to solve (36){(37) in Section 4. We take p = 4 so that these equations can easily be solved by hand. In the course of the calculations, there are two places where square roots must be taken. Thus, there are four possible solutions, depending on whether positive or negative roots are used. If the negative root is always taken, we obtain
dji = (0; : : :; ji; : : :; 0); j = 1; : : :; J; 4
where ji is in the jth position out of J + 1 positions. Set eJi = (0; : : :; 0; Ji). If we let 4
ji
= D dji and 'Ji = D eJi ; 4
4
(13)
and if we assume that D is unitary, then
f j 0; : : :;
j;N=2j ?1g
and f'J 0; : : :; 'J;N=2J ?1g
are orthonormal bases for Wj and BJ , respectively. The union of all these bases is called a wavelet basis for X0 . The signals ji are called wavelets, and the signals 'Ji are called scaling functions. The connection between ji and 'Ji and the DWT output signals is vj (i) = hx0; jii0 and xJ (i) = hx0 ; 'Jii0 : (14) This is easily seen by writing hx0 ; jii0 = hx0 ; D djii0 = (Dx0 ; dji) = (vj ; ji)j = vj (i): The second equation is obtained similarly. Remarks. (i) The Reconstruction Algorithm shows that 1i = H0 1i , and ji
= G0 Gj ?2Hj?1ji; j = 2; : : :; J:
(15)
Also
'Ji = G0 GJ ?1Ji : (16) (ii) Using (11){(12) and (15){(16), it is an easy calculation to show that for j = 1; : : :; J, j;i+l (n)
=
ji (n ? l2j );
p
and 'J;i+l (n) = 'Ji (n ? l2J ); l = 0; : : :; N=2J ? 1: (18) In particular, it follows that the subspace Wj in (9) (resp. BJ in (10)) is closed under circular shifts of 2j (resp. 2J ) units. 3.2. Fast wavelet transforms Let p be an even integer such that 2 p N=2J . As will be shown in Section 4, one can nd a sequence g(0); : : :; g(p ? 1) such that if n = 0; : : :; p ? 1; (19) gj (n) = g(n); 0; n = p; : : :; N=2j ? 1; and if 4
hj (n) = (?1)n gj (1 ? n); n = 0; : : :; N=2j ? 1; 4
p
g(0) = 1 +p 3 ; g(1) = 3 +p 3 ; 4 2 4 2 (21) p p 1 ? 3 ? 3 3 p ; g(3) = p : g(2) = 4 2 4 2 These are the rst four entries in [9, Table 1, p. 980] and [10, Table 6.1, p. 195], where they are called h(n). If we had always taken the positive square root, we would have obtained g~(n) = g(3 ? n), n = 0; 1; 2; 3. The two remaining solutions are ?g(n) and ?g~(n). Let g be given by (21) and de ne gj and hj by (19) and (20) with N = 64 and J = 4. Note that since gj and hj are real, if x0 is real, so is Dx0 . The signals 1;2 , 2;2, 3;2, 4;2, and '4;1 are shown in Figs. 1{5, respectively. The remaining basis signals can be obtained by circular shifts as indicated in (17) and (18). 4. Construction of Gj and Hj for fast wavelet trans-
l = 0; : : :; N=2j ? 1; (17)
forms
(20)
As we saw in Section 3.2, the key to fast wavelet transwhere j = 0; : : :; J ? 1, then the operators Gj and Hj given forms is nding circular FIR lters gj . Then by (20), hj by (11) and (12), respectively, satisfy (4){(7). The discrete will also be circular FIR. Of course, we also need to show 4
that (4){(7) hold for the corresponding operators Gj and Hj . The plan is to construct Gj from G0 and then to construct Hj from Gj . 4.1. Obtaining Hj from Gj Without loss of generality, we may work with G0 and H0 . To simplify the notation, we drop the subscript 0, and we let M = N=2. Then (11) and (12) become
all u. For each u, let g^o(u) be any complex number such that jg^o(u)j2 = 1 ?jg^e(u)j2. Then (27) holds for all u. If we let h(n) be given by (32), then G and H will satisfy (4){ (7). In the next subsection, we consider the more delicate task of constructing g(n) so that in addition to (27), g is also circular FIR. 4.2. Selection of G0 As before we drop the subscript 0. The purpose of NX ?1 this subsection is to nd a circular FIR solution of (5), or g(n ? 2m)x(n); (22) equivalently, (Gx)(m) = (27). n=0 Using (22), it is easy to see that (5) is equivalent to NX ?1 h(n ? 2m)x(n); (23) (Hx)(m) = NX ?1 n=0 g(n + 2m)g(n) = (m); (33) n =0 where g, h, and x are N-periodic sequences. The purpose of this subsection is to use (22){(23) to express (4){(7) in where is M-periodic. Let (m) denote the left-hand side terms of the lters g and h. Proofs of the following lemma of (33). Clearly, has period M. Also, (M ? m) = (m). and theorem, which are adapted from [9, Section 3.A], are Thus, (m) = 0 for m = 1; : : :; M=2 ? 1 if and only if given in Appendix A. (m) = 0 for m = M=2 + 1; : : :; M ? 1. Therefore, (33) holds if and only if (m) = (m) for m = 0; : : :; M=2. We Lemma 1. De ne the M-periodic sequences ge(k) = now impose the constraint g(n) = 0 for p n N ? 1, g(2k); go (k) = g(2k + 1); he(k) = h(2k), and ho(k) = where p M and p is even. It then follows that (33) can h(2k + 1). Denote their M-length DFTs by g^e(u), g^o(u), be replaced by h^ e (u), and ^ho(u), respectively. If G and H are given by pX ?1 (22){(23), then (4) is equivalent to g(n + 2m)g (n) = (m); m = 0; : : :; M=2: (34) 2 2 n =0 ^ jg^e(u)j + jhe(u)j = 1; (24) ^ ^ g^e (u)^go(u) + he (u)ho(u) = 0; (25) Now, in (34), n+2m (p?1)+2(M=2) = p?1+M N ?1, since we assumed p M. Since there is no wraparound, 2 2 ^ jg^o(u)j + jho(u)j = 1; (26) the sum in (34) will automatically be zero if 2m p. Thus, we may rewrite (34) as and (5){(7) are, respectively, equivalent to 4
4
4
4
?1 (27) pX g(n + 2m)g (n) = (m); m = 0; : : :; p=2 ? 1: (35) (28) n=0 (29) Finally, observe that n + 2m p ? 1 if and only if n g^e (u)^he (u) + g^o (u)^ho(u) = 0: p ? 2m ? 1, and since m p=2 ? 1, p ? 2m ? 1 1. Hence, Furthermore, (24){(26) and (29) together imply (35) becomes j^he(u)j = jg^o(u)j and j^ho(u)j = jg^e(u)j: (30) p?X 2m?1 g(n + 2m)g(n) = (m); m = 0; : : :; p=2 ? 1: Similarly, (27){(29) and (25) together also imply (30). In n=0 addition, (27) and (30) imply (24), (26), and (28). (36) Obviously, (36) does not depend on the value of N M Theorem 2 (Wavelet construction). Let G and H be p. Thus, suppose p is even and that we have a solution to given by (22){(23). If (27) holds, and if h is de ned via (36). Suppose also that N and J are such that 2 p J N=2 if gj is given by (19), then (33) will hold h^ e (u) = g^o (u) and ^ho(u) = ?g^e (u); (31) for gj .. IfClearly, we now de ne hj by (20), Theorem 2 shows that (4){(7) will hold for Gj and Hj if they are de ned by (11) then (4){(7) hold. Furthermore, (31) is equivalent to and (12), respectively. As shown in Section 3.2, we obtain (32) a fast wavelet transform. h(n) = (?1)n g(1 ? n): We now observe that (36) gives us only p=2 nonlinear Clearly, this theorem gives a simple recipe to construct equations in p unknowns. Therefore, we must impose adwavelets. Solutions of (27) are easy to construct. Let g^e (u) ditional consistent constraints on g(0); : : :; g(p ? 1) in orbe arbitrary except for being scaled so that jg^e(u)j2 1 for der to x a solution. The linear constraints suggested by
jg^e(u)j2 + jg^o(u)j2 = 1; j^he(u)j2 + j^ho(u)j2 = 1;
5
Daubechies [9] are
but they are not orthonormal unless D is unitary. Now de ne another basis for X0 by setting
p?1 dk X n dz k n=0 g(n)z z=?1 = 0; k = 0; : : :; p=2 ? 1: (37) Daubechies' original reason for introducing (37) was related to regularizing the behavior of the continuous-time analog of (16) for large J [9, Section 3.B]. However, (37) can also be interpreted as saying that the transfer function P corresponding to g, namely pn?=01 g(n)e ?jn! , as well as its p=2 ? 1 derivatives with respect to !, are zero at ! = . In other words, this is a low-pass transfer function with a maximally at stopband. Clearly, the N-length DFT of g(n), denoted by g^(u), is obtained by sampling this transfer function at ! = 2u=N. If h(n) is then de ned by (32), it is easy to see that the corresponding transfer function is high-pass. Theorem 3 (Daubechies). The constraints (37) are consistent with (36). Real-valued solutions g(n) exist, some of which are given in [10, Table 6.1], where they are called h(n). The proof of this result, which is e ectively contained in [9, Section 4.B], is sketched in Appendix A. The theorem statement only guarantees existence of a solution to (36){ (37). In fact, the proof is constructive, obtaining solutions by spectral factorization of certain polynomials; however, this is not the only way to obtain a solution in a given instance. One can attempt a direct solution by hand if the system of equations is small (e.g., p = 4 as in Section 3.3) or by numerical techniques such as Newton's method otherwise. Note that if g(n) is real, then so is h(n). Hence, D can also be viewed as a mapping from IRN to IRN .
= D?1 dji and 'eJi = D?1 eJi : 4
ji
e
4
These signals and the signals in (13) form a biorthogonal set in the sense that h ji ; ej i i0 = hD dji; D?1dj i i0 = (dji; DD?1 dj i ) = ii jj , h ji; 'eJi i0 = h'Ji ; eji i0 = 0, and h'Ji ; 'eJi i0 = ii . From (8) and the de nition of dji and eJi in Section 3.1, we can write 0 0
0 0
0
0 0
0
0
J
N=X 2j ?1
j =1
i=0
X
Dx0 =
0
0
0
vj (i)dji +
N=X 2J ?1 i=0
xJ (i)eJi :
Recalling (14), and noting that x0 = D?1 (Dx0 ), we have the following generalization of (1): x0 =
J
N=X 2j ?1
j =1
i=0
X
hx0 ; jii0
ji
e
+
N=X 2J ?1 i=0
hx0; 'Ji i0 'eJi :
e Let Dx 0 denote the output of the Decomposition Algorithm when Gj is replaced by Aj and Hj is replaced by Cj. e Then De = D?1 and De ?1 = D . Since x0 = De ?1 (Dx 0),
x0 =
J
N=X 2j ?1
j =1
i=0
X
hx0 ; ejii0
ji
+
N=X 2J ?1 i=0
hx0; 'eJi i0 'Ji :
Note also that ji = D dji = D (D eji ) = (D D) eji , and so eji = (D D)?1 ji. This further implies eji = e ?1 e . Of course, ' e e D) (De D) ji Ji and 'eJi ji and ji = (D are similarly related. We say that f ji; 'Ji g and f eji; 'eJig 5. Biorthogonal wavelets are dual wavelet bases (cf. [10, Eq. (3.3.16)]). We now summarize the generalizations of Section 4.1. The wavelets shown in Figs. 1{5 are obviously asymFirst we must augment (22) and (23) with (recall M = metric. Since the Decomposition and Reconstruction AlN=2) gorithms correspond to exact reconstruction lter banks in which the analysis and synthesis lters use the same coeMX ?1 cients, this asymmetry is to be expected [22]. The diculty a(n ? 2m)x(m); (Ax)(n) = arises because we were too restrictive in Section 2.1 when m =0 we required T T = I and TT = I. We can relax these MX ?1 conditions by requiring only that T be invertible, where c(n ? 2m)v(m); (Cv)(n) = T ?1 has the form T ?1 = [A C ]. Then the generalizam=0 tions of (4){(7) are AG + CH = I, GA = I, HC = I, and GC = 0, respectively. We also need to add HA = 0. where the complex conjugate is included only for conveThe Decomposition Algorithm, which uses only Gj and nience. After reading the proof of Lemma 1 in Appendix A, Hj , stays the same, and its output is still called Dx0. The it is easy to see that in the present context, AG + CH = I Reconstruction Algorithm, which uses only Gj and Hj , is equivalent to still computes D , but now D 6= D?1 . However, we can (38) g^e(u)^ae (u) + ^he (u)^ce (u) = 1; compute D?1 if we modify the Reconstruction Algorithm by replacing Gj by Aj and Hj by Cj , assuming that Aj Gj + ^ (39) g^e (u)^ao(u) + he (u)^co (u) = 0; Cj Hj = I for each j. ^ g^o(u)^ao (u) + ho (u)^co (u) = 1; (40) If D is invertible, then so is D . In this case, the signals N ^ j g^o (u)^ae(u) + ho (u)^ce (u) = 0; (41) ji and 'Ji de ned in (13) still form a basis for X0 = C , 6
where (p; p~; m) = minfp ? 2m ? 1; p~ ? 1g. (If p~ = 1 or 2, we have only the rst equation.) Clearly, these conditions will hold for all N M maxfp; p~; 2g. 5.1. Symmetric, biorthogonal, FIR wavelets Let ( (maxfp; p~g ? 1)=2; p; p~ both odd; K = maxfp; p~g=2 ? 1; p; p~ both even: We view both g(n) and a(n) as being of length 2K + 1 in the rst case, and of length 2K + 2 in the second case. In other words, we pad the tail of the shorter sequence with zeros. We now require that g and a be symmetric. For p and p~ both odd, this means that g(K + n) = g(K ? n) and a(K + n) = a(K ? n) for n = 1; : : :; K. For p and p~ both even, this means that g(K +n) = g([K +1] ? n) and a(K + n) = a([K+1] ? n) for n = 1; : : :; K+1. Note that since we zero-padded the tail of the shorter sequence, by imposing symmetry, we are forcing zeros in the head as well. To pin down choices further for g and a, we impose conditions analogous to (37). With slight abuse of notation, let
and GA = I, HC = I, GC = 0, HA = 0 are, respectively, equivalent to g^e (u)^ae (u) + g^o(u)^ao (u) = 1; (42) ^he (u)^ce (u) + h^ o (u)^co (u) = 1; (43) g^e (u)^ce (u) + g^o (u)^co (u) = 0; (44) ^he (u)^ae (u) + h^ o (u)^ao (u) = 0: (45) Suppose that g and a are such that (42) holds. Taking our hint from (31), let h and c be de ned via h^ e (u) = a^o (u) and h^ o(u) = ?a^e (u); (46)
4
4
and
c^e (u) = g^o(u) and c^o (u) = ?g^e (u): (47) Then (38){(41) and (43){(45) also hold. Note that (46) and (47) are equivalent to h(n) = (?1)n a(1 ? n) and c(n) = (?1)n g(1 ? n); (48) respectively. Thus, we now have a simple recipe to construct biorthogonal wavelets. If g^e (u) and g^o (u) are arbitrary except for the constraint that g^e (u)^go(u) 6= 0 for all u, then we can always nd choices of a^e (u) and a^o(u) that satisfy (42). If h and c are then de ned by (48), and if we can do this for j = 0; : : :; J ? 1, then we have a system of biorthogonal wavelets. In order that the algorithms for computing D, D , and ? 1 D be fast, we would like to construct compactly supported biorthogonal wavelets based on two sequences, to be chosen below, say g(0); : : :; g(p ? 1) and a(0); : : :; a(~p ? 1), where we assume that 1 minfp; p~g and maxfp; p~; 2g N=2J . We do not require that p or p~ be even. We de ne the N=2j -periodic sequence gj by (19), and de ne aj similarly. We use (48) to de ne hj in terms of aj and cj in terms of gj . To be sure this will all work correctly, we must generalize the derivation leading from (33) to (36). Now, GA = I and (42) are both equivalent to NX ?1
G(z) = 4
pX ?1 n=0
g(n)z n and A(z) = 4
pX ~?1 n=0
a(n)z n : (51)
We require that the derivatives with respect to z satisfy G(k)(?1) = 0; k = 0; : : :; ` ? 1; (52) A(k) (?1) = 0; k = 0; : : :; `~ ? 1; where 1 ` p and 1 `~ p~. Note that these conditions do not uniquely specify a solution of (49); if g(n) and a(n) satisfy (50) and (52), then so do g(n)= and a(n) for any nonzero constant . We therefore impose the normalization condition p~?1 X p (53) a(n) = 2:
p
n=0
The choice of 2 is only for convenience; any nonzero con(49) stant would do. n=0 Example. Let p = 6 and p~ = 4. Then K = 2. By imposwhere g and a are regarded as N-periodic and is M = ing symmetry, only a(3) = a(2), g(5) = g(0), g(4) = g(1), N=2-periodic. Let (m) denote the left-hand side of (49). and g(3) = g(2) are nonzero. Substituting these into (50) Observe that (m) = (m) for m = 0; : : :; M ? 1 if and only yields the two independent equations, g(2)a(2) = 12 and if (m) = (m) for m = 0; : : :; M=2 and (M ? m) = 0 g(0) = ?g(1). Take ` = 2 and put no constraint on A(z). for m = 1; : : :; M=2 ? 1. If we now impose the conditions Since g is symmetric, G(0)(?1) = 0 already. G0 (?1) = 0 g(n) = 0 for n = p; : : :; N ? 1 and a(n) = 0 for n = implies p g(2) = 8g(1). The normalization (53) implies p~; : : :; N ? 1, and proceed as in Section 4.2, we see that a(2) = p2=2. Thus, g(2) = p2=2, g(1) = p2=16, and (49) is equivalent to the pair of equations, g(0) = ? 2=16. The remaining coecients are pfound from symmetry. If we divide all our coecients by 2, we ob(p; p~;m) X tain the coecients in the second row of Table 8.2 in [10, g(n + 2m)a(n) = (m); 0 m < p=2; p. 277]. For N = 64 and J = 3, plots of e3;2, 'e3;1, 3;2, n=0 (50) and '3;1 are shown in Figs. 6{9, respectively. (~pX ;p;m) When we imposed (52), there was no guarantee that g(n)a(n + 2m) = 0; 1 m < p~=2; these constraints were consistent with (42). That this is n=0 g(n + 2m)a(n) = (m);
7
indeed the case essentially follows from [10, Section 8.3.4]. Recall that (36) was derived from (33), and (33) is We brie y sketch the argument in Appendix A. equivalent to (27). Letting g^(u) denote the N-length DFT of g(n), and using the analog of (55), we nd that (27) is Appendix A equivalent to Proof of Lemma 1. To show that (4) is equivalent to jg^(u)j2 + jg^(u + N=2)j2 = 2: (57) (24){(26), rst note that (22){(23) imply that (4) is equivalent to Now observe that g^(u) = G(e j! ) when ! = ?2u=N and u is an integer. If (56) holds, then MX ?1 g(n ? 2m)g( ? 2m)+ h (n ? 2m)h( ? 2m) = (n ? ); jG(e j! )j2 = [cos2(!=2)]p=2 jF(e j! )j2: m=0 (54) Next, if F(z) has real coecients, then jF(e j! )j2 also has where denotes the N-periodic Kronecker delta. By con- real coecients and is real-valued. Hence, we can write sidering (54) when n and are both even, it is easy to jF(e j! )j2 = F(e j! )F(e ?j! ) as a polynomial of degree P ?1 see that (54) implies mM=0 ge(k ? m)ge (?m) + he(k ? p=2 ? 1 in cos !, or equivalently, in sin2 (!=2); i.e., we may m)he (?m) = (k), where here is M-periodic. Tak- write jF(e j! )j2 = P(sin2 (!=2)), where P is a polynomial ing an M-length DFT and changing ?u to u yields (24). of degree p=2 ? 1. So, if (56) holds for some F with real Similarly, taking n odd and even in (54) leads to (25). coecients, and if (57) holds (which is what we are trying Taking n and both odd leads to (26). To show that to prove), then (5){(7) are, respectively, equivalent to (27){(29), proceed y p=2P(1 ? y) + (1 ? y)p=2 P(y) = 2; (58) as follows. We treat only (5) as the other two are similar. As noted in Section 4.2, (5) is equivalent to (33). where y = cos2 (!=2), ! = ?2u=N, and u is an integer. Breaking the sum in (33) into even and odd n yields PM ?1 Conversely, the polynomial k=0 ge(k+m)ge (k)+go (k+m)go (k) = (m). Taking an M-length DFT of this equation gives (27). The remaining p=X 2?1 p=2 ? 1 + k statements of the theorem are easily veri ed. 2 yk (59) P(y) = 2 k k =0 Proof of Theorem 2. First note that (31) trivially implies (25) and (29). Next, (31) implies (30), and by hy- satis es (58).1 Furthermore, P(y) is nonnegative for 0 pothesis, (27) holds. Hence, (24), (26), and (28) hold by y 1, and by Riesz's lemma, there exists a (nonunique) Lemma 1. Since (24){(29) now hold, Lemma 1 further im- polynomial Q of degree p=2 ? 1 with real coecients that plies that (4){(7) hold. To show that (31) and (32) are satis es jQ(e j! )j2 = P(sin2(!=2)). Thus, if we take the Q equivalent, substitute (31) into the easily veri ed formula, of Daubechies2 and obtain the g(n) from ^h(u) = ^he (u) + h^ o (u)e ?j2u=N ; G(z) = [ 21 (1 + z)]p=2 Q(z); (55) where h^ is the N-length DFT of h. We obtain h^ (u) = and set ! = ?2u=N and y = cos2 (!=2), then g^o (u) ? g^e (u)e ?j2u=N , from which (32) follows. The fact jg^(u)j2 + jg^(u + N=2)j2 that (32) implies (31) follows by separately considering (32) = jG(e j! )j2 + jG(e j(!?) )j2 with n replaced by 2n + 1 and by 2n. 2 = y p=2P(1 ? y) + (1 ? y)p=2 P(y) Proof of Theorem 3. The proof is self-contained ex= 2: 2 cept for the reference to Riesz's lemma (see e.g., [9, p. 972, Lemma 4.2] or [10, Lemma 6.1.3]), which is a result about Proof of Consistency of (52). We treat the case in the spectral factorization of polynomials that can be read which p and p~ are both even, the other case being similar. independently of the rest of [9] or [10]. First note that just as (57) was established above, (42) is Let G(z) be as in (51). Since (37) is equivalent to saying equivalent to that (1 + z)p=2 is a factor of G(z), g^(u)^a(u) + g^(u + N=2)^a(u + N=2) = 2: (56) G(z) = [ 21 (1 + z)]p=2F(z) 1 ?p=2 p=2 4
4
4
If P (y) satis es (58), P (y) = (1 ? y) [2 ? y P (1 ? y)]. Substituting the Taylor series of (1 ? y)?p=2 yields (59) plus higherorder terms. Since P (y) is a polynomial of degree less than p=2, the higher-order terms must sum to zero. Eq. (59) is the unique lowestdegree solution of (58). Higher-degree solutions also exist [10, p. 171].
for some polynomial F of degree p=2 ? 1. The factor of 1=2 is included for convenience in the derivation below. Observe that the coecients of G are all real if and only if the coecients of F are all real. We now require that 2 In general there are p=2 choices for Q(z ) (and hence G(z )). The g(n) be real. Thus, it suces to prove the existence of a coecients (n) (called h(n) in [9], [10]) given in [9, p. 980, Table 1], polynomial F of degree p=2 ? 1 with real coecients such [10, p. 195, gTable 6.1] correspond to taking Q(z) to have all its roots that if G(z) is de ned by (56), then (36) holds. inside the unit circle; i.e., Q(z) is a minimum-phase lter. 8
If we set ! = ?2u=N, this can be written as G(e j! )A(e j! ) + G(e j(!?) )A(e j(!?) ) = 2:
[5] G. Caire, R. L. Grossman and H. V. Poor, \Wavelet transforms associated with nite cyclic groups", IEEE Trans. Inform. Theory, Vol. 39, No. 4, July 1993, pp. 1157{1166. [6] C. K. Chui, An Introduction to Wavelets, Academic Press, New York, 1992. [7] A. Cohen, I. Daubechies and J.-C. Feauveau, \Biorthogonal bases of compactly supported wavelets", Commun. Pure Appl. Math., Vol. 45, 1992, pp. 485{560. [8] R. R. Coifman and M. V. Wickerhauser, \Best-adapted wavelet packet bases", preprint, February 1990. [9] I. Daubechies, \Orthonormal bases of compactly supported wavelets", Commun. Pure Appl. Math., Vol. 41, November 1988, pp. 909{996. [10] I. Daubechies, Ten Lectures on Wavelets, SIAM, Philadelphia, PA, 1992. [11] I. Gohberg and S. Goldberg, Basic Operator Theory, Birkhauser, Boston, 1980. [12] C. E. Heil and D. F. Walnut, \Continuous and discrete wavelet transforms", SIAM Rev., Vol. 31, No. 4, December 1989, pp. 628{666. [13] F. Hlawatsch and G. F. Boudreaux-Bartels, \Linear and quadratic time-frequency signal representations", IEEE Signal Process. Mag., Vol. 9, No. 27, April 1992, pp. 21{67. [14] K. Ho man and R. Kunze, Linear Algebra, 2nd Edition, Prentice-Hall, Englewood Cli s, NJ, 1971. [15] S. Mallat, \Multiresolution approximations and wavelet orthonormal bases of L2 (IR)", Trans. Amer. Math. Soc., Vol. 135, No. 1, September 1989, pp. 69{88. [16] Z. J. Mou and P. Duhamel, \Short-length FIR lters and their use in fast nonrecursive ltering", IEEE Trans. Signal Process., Vol. 39, June 1991, pp. 1322{1332. [17] W. H. Press, S. A. Teukolsky, W. T. Vetterling and B. P. Flannery, Numerical Recipes in C: The Art of Scienti c Computing, 2nd Edition, Cambridge Univ. Press, New York, 1992. [18] O. Rioul, \A discrete-time multiresolution theory", IEEE Trans. Signal Process., Vol. 41, No. 8, August 1993, pp. 2591{ 2606. [19] O. Rioul and P. Duhamel, \Fast wavelet algorithms for discrete and continuouswavelet transforms", IEEE Trans. Inform. Theory, Vol. 38, No. 2, March 1992, pp. 569{586. [20] O. Rioul and M. Vetterli, \Wavelets and signal processing", IEEE Signal Process. Mag., Vol. 8, No. 4, October 1991, pp. 14{ 38. [21] M. J. Shensa, \The discrete wavelet transform: Wedding the a trous and Mallat algorithms", IEEE Trans. Signal Process., Vol. 40, No. 10, October 1992, pp. 2464{2482. [22] M. J. T. Smith and T. P. Barnwell III, \Exact reconstruction techniques for tree-structured subband coders", IEEE Trans. Acoust. Speech, Signal Process., Vol. ASSP-34, No. 3, June 1986, pp. 434{441. [23] G. Strang, \Wavelets and dilation equations: A brief introduction", SIAM Rev., Vol. 31, No. 4, December 1989, pp. 614{627. [24] G. Strang, \Wavelet transforms versus Fourier transforms", Bull. Amer. Math. Soc., Vol. 28, No. 2, April 1993, pp. 288{ 305. [25] M. Vetterli, \Running FIR and IIR ltering using multirate lter banks", IEEE Trans. Acoust. Speech, Signal Process., Vol. ASSP-36, May 1988, pp. 730{738. [26] M. Vetterli and C. Herley, \Wavelets and lter banks: Theory and design", IEEE Trans. Signal Process., Vol. 40, No. 9, September 1992, pp. 2207{2232.
(60)
Now, the symmetry of g(n) and a(n) implies G(e j! ) = e j!(K +1=2) cos(!=2) q(cos !); A(e j! ) = e j!(K +1=2) cos(!=2) r(cos !);
(61)
where q and r are polynomials of degree less than or equal to K. If (52) also holds, then G(z) = (1 + z)` G0(z) and A(z) = (1 + z)`~A0(z) (62) for some polynomials G0 (z) and A0 (z). It can be shown that (61){(62) imply that G(e j! ) = e j!(K +1=2) cos+1 (!=2) q0(cos !); A(e j! ) = e j!(K +1=2) cos~+1 (!=2) r0 (cos !); where q0 and r0 are polynomials, and where = ` if ` is even and = ` ? 1 if ` is odd and similarly for ~. Substituting this into (60), we obtain cos2L(!=2) q0 (cos !) r0 (cos !) + sin2L (!=2) q0(? cos !) r0(? cos !) = 2; where L = ( + ~ + 2)=2 is a positive integer. Since cos ! = 1 ? 2 sin2(!=2), we can write q0(cos !)r0(cos !) = P(sin2 (!=2)) for some polynomial P. Letting y = sin2 (!=2), we have (1 ? y)L P(y) + yL P(1 ? y) = 2: As shown above in the proof of Theorem 3, there are polynomials P that satisfy this equation. By factoring P, feasible q0 and r0 can be found. Note that the choice of factors imposes a normalization condition, though not necessarily the one in (53). 2
Acknowledgements
This work was supported in part by the Air Force Oce of Scienti c Research under Grants AFOSR{90{0181 and F49620{92{J{0305.
References
[1] B. Alpert, G. Beylkin, R. Coifman and V. Rokhlin, \Waveletlike bases for the fast solution of second-kind integral equations", SIAM J. Sci. Comput., Vol. 14, No. 1, January 1993, pp. 159{184. [2] G. Beylkin, \On the representation of operators in bases of compactly supported wavelets", preprint. [3] G. Beylkin, R. Coifman and V. Rokhlin, \Fast wavelet transforms and numerical algorithms I", Comm. Pure Appl. Math., Vol. 44, 1991, pp. 141{183. [4] P. J. Burt and E. H. Adelson, \The Laplacian Pyramid as a compact image code", IEEE Trans. Commun., Vol. COM-31, April 1983, pp. 532{540.
9
1 0.9 0.8 0.7 0.6 0.5 0.4
0.5
0.3 0.4
0.2 0.1
0.3
0 -0.1
0.2
-0.2 0.1
-0.3 -0.4
0
-0.5 0
10
20
30
40
50
-0.1
60
-0.2
Fig. 1. Wavelet 1;2(n) of Section 3.3.
-0.3 -0.4 0
10
20
30
40
50
60
Fig. 4. Wavelet 4;2(n) of Section 3.3.
0.8 0.7 0.6 0.5 0.4 0.3 0.2 0.1 0 -0.1 -0.2 -0.3 -0.4 -0.5 -0.6 0
10
20
30
40
50
60
Fig. 2. Wavelet 2;2(n) of Section 3.3.
0.40 0.35 0.30 0.25
0.6
0.20
0.5
0.15
0.4
0.10
0.3 0.05 0.2 0.1
0
0
-0.05
-0.1
-0.10
-0.2
0
10
20
30
40
50
60
-0.3
Fig. 5. Scaling function '4;1(n) of Section 3.3.
-0.4 -0.5 -0.6 0
10
20
30
40
50
60
Fig. 3. Wavelet 3;2(n) of Section 3.3. 10
0.4 0.3 0.2 0.1 0 -0.1 -0.2 -0.3 -0.4 0
10
20
30
40
50
60
Fig. 6. Biorthogonal wavelet e3;2(n) of Section 5.
0.45
0.40
0.40
0.35
0.35 0.30
0.30
0.25
0.25 0.20
0.20 0.15 0.15
0.10
0.10
0.05 0
0.05
-0.05
0
-0.10 0
10
20
30
40
50
0
60
10
20
30
40
50
60
Fig. 7. Biorthogonal scaling function 'e3;1(n) of Section 5. Fig. 9. Biorthogonal scaling function '3;1(n) of Section 5.
0.5 0.4 0.3 0.2 0.1 0 -0.1 -0.2 -0.3 -0.4 -0.5 0
10
20
30
40
50
60
Fig. 8. Biorthogonal wavelet 3;2(n) of Section 5. 11
