mathematics of computation volume 63,number 208 october 1994,pages 689-699 ON THE COMPUTATION OF BATTLE-LEMARIE'S WAVELETS MING-JUN LAI Abstract. We propose a matrix approach to the computation of BattleLemarié's wavelets. The Fourier transform of the scaling function is the product of the inverse F(x) of a square root of a positive trigonometric polynomial and the Fourier transform of a B-spline of order m . The polynomial is the symbol of a bi-infinite matrix B associated with a B-spline of order 2m . We approximate this bi-infinite matrix B2m by its finite section As , a square matrix of finite order. We use As to compute an approximation \s of x whose discrete Fourier transform is F(x). We show that xs converges pointwise to x exponentially fast. This gives a feasible method to compute the scaling function for any given tolerance. Similarly, this method can be used to compute the wavelets. 1. Introduction Battle-Lemarié's wavelets [1,3] may be constructed by using a multiresolution approximation built from polynomial splines of order m > 0. See, e.g., [4] or [2]. To be precise, let Vq be the vector space of all functions of L2(R) which are m - 2 times continuously differentiable and equal to a polynomial of degree m - 1 on each interval [« + m/2, n+l + m/2] for all « € Z. Define the other resolution space Vk by Vk:={u{2kt): u£ V0}, VfceZ. It is known that {Vk}k£z provide a multiresolution approximation, and there exists a unique scaling function <p such that Vk= %oanLi{2kl2fa2kt- ri) : n € Z} for all k , and the integer translates of tp are orthonormal to each other. (See, e-g->[4]-) Define a transfer function //(w) by ffim\ = - —;—r-, ^(2ft)) H{0}) <pi(0) where tp denotes the Fourier transform of <p. Then the wavelet y/ associated with the scaling function tp is given in terms of its Fourier transform by Vita) = e-j0}/2Hico/2+ n)faco/2). Received by the editor May 14, 1993. 1991MathematicsSubjectClassification.Primary 41A15,41A30,42C15, 47B35. Key words and phrases. B-spline, bi-infinite matrices, exponential decay, finite section, positive operator, Toeplitz matrix, wavelet. Research supported by the National Science Foundation under Grant #DMS-9 303121. ©1994 American Mathematical Society 0025-5718/94$1.00+ $.25 perpage 689 License or copyright restrictions may apply to redistribution; see https://www.ams.org/journal-terms-of-use 690 MING-JUN LAI Here and throughout, ; := v^-T • The scaling function tp associated with the multiresolution approximation may be given by (1) faca) = ^ l Bmico), \lïlkzz\Bm(<o+ 2kn)\* where Bm is the well-known central B-spline of order m whose Fourier trans- form is given by g . . fsinco/2' By using Poisson's summation formula, we have faœ) = V/£*ez*2m(fc)e-*0' Bmiœ) Thus, the transfer function is <2) g""-fet^(~B/2)-- Then the wavelet y/ associated with <p is given by (3) ytiœ) = e-J»'2Hiœ/2 + n) l ^EkezB2m(k)e-^/2 Bmioo/2). The above Fourier transforms of <p, H, and y/ suggest that the scaling function, transfer function, and wavelet have the following representations: <PÍt)= YakBm(t-k), kez Hi(o)= Y,ßke-ik,°, kez v(t) = YykBm(2t-k). kÇZ In this paper, we propose a matrix method to compute the a^'s, ßk's, and yks. Let us use <pto illustrate our method as follows: view "51ke7jB2mik)e~ikb} as the symbol of a bi-infinite matrix Ti2m= ibik)itkez with b¡ k = ¿o,fc-i = B2m{k-i) for all i,k£Z. Similarly, JY,k€ZB2m(k)e-Jk(0 can be viewed as the symbol of another (unknown) bi-infinite matrix C2m. Then it is easy to see that To find kez is equivalent to solving JZkezB2m(k)e-Jk» C2mx = S with S = (<5,),ez, So = I, and S¡ = 0 for all i £ Z\{0}, where x = (a^)fceZ. Our numerical method is to find an approximation to x within a given tolerance. License or copyright restrictions may apply to redistribution; see https://www.ams.org/journal-terms-of-use ON THE COMPUTATIONOF BATTLE-LEMARIÉ'S WAVELETS 691 Let An = ib¡k)-N<i,k<N oe a finite section of B2m. Note that AN is symmetric and totally positive. Thus, we can find P# such that H = AN by using, e.g., the singular value decomposition. Then we solve ?n*n = Sn with <5jva vector of 2N+1 components which are all zeros except for the middle one, which is 1. We can show that x^ converges pointwise to x exponentially fast. Similarly, we can use this idea to compute an approximation of {ßk}kcz by (2) and {Yk}kez hy (3). Therefore, the discussion mentioned above furnishes a numerical method to compute Battle-Lemarié's wavelet. To prove the convergence of xjv to x, we place ourselves in a more general setting. We study a general bi-infinite matrix A . (For the case of BattleLemarié's wavelets, A = B2m.) We look for certain conditions on A such that the solution x^ of PnXn = Sn with PN = An converges to the solution x of Px = S with P2 = A, where An is a finite section of A. This is discussed in the next section. In the last section, we show that the bi-infinite matrix H2m satisfies the conditions on A obtained in §2. This will establish our numerical method for computing Battle-Lemarié's wavelets. 2. Main results Let Z be the set of all integers. Let I2 := 12{L) be the space of all square summable sequences with indices in Z. That is, /2(Z)= U... ,x-x,Xq,Xx, ...)' : Y M2<°°L I «=—00 ) It is known that I2 is a Hubert space. We shall use x to denote each vector in I2 and use A to denote a linear operator from I2 to I2. It is known that A can be expressed as a bi-infinite matrix. Thus, we shall write A = iaik)ikez ■ In the following, we shall consider A to be a banded and/or Toeplitz matrix. That is, A is said to be banded if there exists a positive integer b such that aik = 0 whenever \i—k\ > b. The matrix A is said to be Toeplitz if ai+ktm+k = a^m for all i, k, m £ Z. Denote by F{x.)ioj) the symbol of a vector x £ I2, i.e., Fix)iœ) = 1£x'e~jiœ¡ez Denote by FiA)ico) the symbol of a Toeplitz matrix A = iaik)itk€Z >i.e., FiA)i(o) = Ydai,oe-»°>. /ez Suppose that FiA)ia) ¿ 0 and £;€zla',o| < oo. It is known from the well-known Wiener's theorem that there exists a sequence x such that with Yjk \xk\ < °° • It is easy to see that to find this sequence x is equivalent to solving the linear system of bi-infinite order: Ax = ô, License or copyright restrictions may apply to redistribution; see https://www.ams.org/journal-terms-of-use MING-JUN LAI 692 where Ô = (..., ¿_i, S0, Sx, ... )' with ô0 = 1 and <5,= 0 for all î e Z\{0} . Furthermore, if the matrix A is a positive operator, then there exists a unique positive square root P of A . That is, P2 = A . The symbol representation is FiP)ico) = y/FiA)ico). To find /r(P)(w) is equivalent to finding a matrix P such that P2 = A . Certainly, we cannot solve a linear system of bi-infinite order. Neither can we decompose a matrix of bi-infinite order into two matrices of bi-infinite order. However, we can do this approximatively. Let TVbe a positive integer, and let An = (Qik)-N<i,k<N be a finite section of A. Let In,oo = (0, /2jv+i,2jv+i . 0) be a matrix of 2TV+ 1 rows and bi-infinite columns with I2n+x,2n+\ being the identity matrix of size (27V+ 1) x (27V+ 1) such that An = 'N ,ooAIn ,oo • Denote Sn = In,ooS and x^ = IN>00x. Then we shall solve the following linear system: AnXn = Sn ■ We claim that xN converges to x exponentially fast as TV increases to oo, under certain conditions on A . Furthermore, we shall solve PN = An for PN by using the singular value decomposition. Once we have Pn , we shall solve PnSn —Sn ■ We claim that y^ converges to y exponentially fast as TV—»oo, provided A satisfies certain conditions. To check the conditions on A, we need the following definition. Definition. A matrix A = (a,fc)i,itez is said to be of exponential decay off its diagonal if \aik\ < K^-k\ for some constant K and r £ (0, 1). We begin with an elementary lemma. Lemma 1. Suppose that A is of exponential decay off its diagonal and has a bounded inverse. Suppose that ANX= iaik)-N<i,k<N satisfies the property that |âJjfc(7V)|<A>l'-*l, V-N<i,k<N, for all TV> 0. Then there exists rx £ (0, 1) and a constant Kx such that Wlff^x-XNh^Kxrf, where x is the solution of Ax = 5 and xN is the solution of AnXn = Sn ■ Proof. From the assumption of the lemma, there exist K and re(0, 1) such that A = (a/*),,fcez and ANl = (àiik(N))_N<i>ksN satisfy lOfltl< Jfr"-*l and \âik\< Kr^ . Write B AIn,oo — AN C and d = BANXSN with d = {.. , d-N-x, d-N)'• License or copyright restrictions may apply to redistribution; see https://www.ams.org/journal-terms-of-use ON THE COMPUTATIONOF BATTLE-LEMARIE'S WAVELETS 693 Then we have, for each i = -co, ... , -TV - 1, -TV, N \di\= Y aikàk,o(N) N <K2 £ rXi-k\r\k\ \k=-N ■K2 ír-'¿r2/c + TVr-'J < CX~l for some constant C and A £ (0, 1). Thus, ||iL4_1<5jv||2 < C'kN. Similarly, \\CANXÔN\\2<C'XN. Hence, ||/v,oox - xN\\2 < ||x - IN>0oxN\\2 < \\A-X\\2\\Ô - AIn¡00AnxOn\\2 ' B ' <U- s- An AN SN C 2 BA-NX <\\A- hN+X,2N+l CANX Sn < \\A-x\\2í\\BA-xOn\\2+ \\CA^xSn\\2) < ||^-1||22C'A7V, hence the assertion with Kx = 2C||^_1||2 lemma. D and rx = X. This establishes the Next, we consider approximating the square root of a positive operator. Lemma 2. Let P be the unique square root of a positive operator A. Suppose that A is banded and \A - I\2 < r < 1, where I is the identity operator from I2 to I2. Then P = ÍPik)i,kez is of exponential decay off its diagonal. Proof. The uniqueness of P and the convergence of the series ,(2_/-3)!! i=0 (2/)!! {A l) imply that P = ^=^I + iA-i) = ^i-iyVl J>J\a-i) ;=0 *• '" The matrix A is banded and so is A - I. If A - I has bandwidth b, then iA - I)1 is also banded with bandwidth ib . Thus, \pik\ < Kr^'~k^b for some constant K. This finishes the proof. D Lemma 3. Let P be the unique square root of a positive operator A. Suppose that A is banded and \\A - I\\2 < r < I, where I is the identity operator from I2 to I2. Then P~x = ÍPik)¡,kez is of exponential decay off its diagonal. Proof. The uniqueness of P~x and the convergence of the series £(-D /=0 t(2f-l)ll (20« (A-iy License or copyright restrictions may apply to redistribution; see https://www.ams.org/journal-terms-of-use 694 MING-JUN LAI imply that ¿(2/-1)!! /»-i = iA)~xl2 = (I + (A- I))'1'2 = £(-1) U-/y (2/)!! ¡=o Now we use the same argument as in the lemma above to conclude that P~x is of exponential decay off its diagonal. D Let Pff be the square root of AN. That is, PN = An ■ Denote Pn = În,ooPIn,oo • ^e neec* t° estimate PnPn - PnPn ■ We have Lemma 4. Let R = (rik)-N<iik<N := PnPn - PnPn ■ Then rik = 0(rAr/(4*))for k = -TV/4+ 1, ... , TV/4- 1 and i = -TV, ... , TV,where b is the bandwidth of A and r is as defined in Lemma 3. Proof. It is known that P and A commute. Let us write ax P= B' B a2 C C*4 PN C «3 and A= 'ßi a a1 An ß2 c' fo c &. We have B'a + PNAN + C'c = a'B + ANPN + c'C. Thus, PNAN - ANPN = a'B - B'a + c'C - C'c. Let E = a'B - B'a + c'C - C'c and IN := 72jv+i,2n+x • We have PnÍAn - In) = ÍAn - In)Pn + E and n-X PnÍAn - IN)n = (AN- InYPn + £(4v - IN)kEiAN - In)"^1 k=0 by using induction. Then, we have PnPn „(2*-3)11 PnÍAn - In)" n=0 v ' _yv 1)H(2«-3)!! (2«)!! ÍAn - In)"Pn «=o oo ,~ + Y^-^-^w1 „=o "5\|l n—X £(^ ^n)" PnPn + £(-l)"k k=0 ,„(2n -3)!! ,' _n \¿n)~ - InYeíAn - iN)n-k-x n-X /" £(^v ._„ - IN)kEiAN - IN)n- A;-l To estimate T? = P#^V - PnPn which is the summation above, we break R into two parts and estimate the first by 1MI n-X O (2«) /i=Vl + l oo * Y A:-l fc=0 n- (2"(2«)!!! " IIEWAn~ In11"2 - *' IIAn~ In^ J2n) n=NX+X ^ '" Thus, this part has the desired property if we choose TV1appropriately. Next, we note that An-In is banded and its bandwidth is b. Thus, for 0 < « < TV1, ÍAn - In)" is also banded and has bandwidth nb < bNl. License or copyright restrictions may apply to redistribution; see https://www.ams.org/journal-terms-of-use ON THE COMPUTATIONOF BATTLE-LEMARIÉ'S WAVELETS 695 Note also E = ieik)_N<i,k<N has the following property: JO for -N + b<k<N-b,-N eik ~ 1 0(^-1*1) for -N<i<-N + b<i<N-b, + b and N-b<i<N, -N<k<N. It follows that iAN - In)1E has a similar property as E : 0 for -N + b<k<N-b, -N + kb + b<i<N-lb-b, ((AN - IN)'E)ik = < 0(r*-l*l) for -N<i<-N + lb + b, N-lb-b<i<N, and -N<k<N. Choose TV1such that TV/(4Z>) < TV1< TV/(4¿>) + 1. Then iAN - I)NX has bandwidth ¿>TV1 < TV/4+ b and hence «A„-r)<E(A„-r)"-'-% = ( °('3""4-'-W) I 0(1) if 1*1S "/* and - A-< i < AT, otherwise for / = 1, ... , TV1. Putting these two parts together, we have established that R has the desired property. D We are now ready to prove the following. Theorem 1. Suppose that A is a positive operator and \\A - I\\2 < 1. Suppose that A is a banded matrix. Let P be the unique square root of A and y the solution of Py = ô. Let PN be a square root matrix such that PN = AN and $n the solution of PnSn —Sn ■ Then ||Tv,ooy-yv||2 <KXN for some X £ (0, 1) and a constant K > 0. Proof. Let P = ip¡k)i,k£z and Pn = ÍP¡k)-N<i,k<N■ By Lemma 2, the matrix P is of exponential decay off its diagonal. By Lemma 3, we know that Pn is of exponential decay off its diagonal uniformly with respect to TVbecause of \\An-I2n+x,2n+xII2 < 1, which follows from \\A-1\\2 < 1. The invertibility of A implies that P is invertible. From ||v4-/||2 < 1 it follows that the inverse of P is bounded. Let y^ be the solution of PnVn = Sn ■ Thus, we apply Lemma 1 to conclude that ll/jv.ooy-yvlb <KxrN for some r e (0, 1). We now proceed to estimate ||yiv - y^lb • Note that P2 = A implies AN = PN + B'B + C'C or P2-P2= Thus, we have iPN + PN)iPN - Pn) = Pn-Pn B'B + CC. + PnPn - PnPn = B'B + CC + R, where R was defined in Lemma 4. Hence, iPN - PN) = iPN + Pn)-\B'B + C'C + R). Note that the entries of B'B + C'C have the exponential decay property: iB'B + C'C)ik = 0{rN~\k\). By Lemma 4, we know that each entry of the middle section (TV/2 columns) of the columns of B'B + C'C + R has exponential License or copyright restrictions may apply to redistribution; see https://www.ams.org/journal-terms-of-use 696 MING-JUN LAI decay 0(rAr/<4ft>). Both PN and PN are positive and ||(¿V + Pa,)-1||2 < ll^'lb is bounded. Recall that P^x is of exponential decay off its diagonal. We have llyv-yv||2 < II^IHI^ - PnPñ1Sn\\2 <\\PÑx\\2\\(Pn-Pn)(PñxSn)\\2 < \\PñX\\2\\(Pn+ Pn)-1\\2\\(B'B + C'C + R)P^xSn\\2 <KXN for some X £ (r, 1). This completes the proof. D In the proof above, an essential step is to show that each entry of the middle section of the columns of Pn - Pn is of exponential decay. This indeed follows from i?N - PN) = (Pn + Pn)~x(B'B + C'C + R), the boundedness of (Pn + TV)-1 , and the fact that each entry of the middle section of the columns of B'B + CC + R is of exponential decay. This has its own interest. Thus, we have the following Theorem 2. Suppose that A is a positive operator and \\A - I\\2 < 1. Suppose that A is a banded matrix. Let P be the unique square root of A and Pn = In,ooP(In,ooY • Let Pn be a square root matrix such that PN = AN ■ Then \\PNSN- PnSn\\2 < KXN for some X £ (0, 1) and a constant K. Finally, we remark that if ||^-/||2= 1, then each entry of the middle section of the columns of R is convergent to 0 with speed ^ . The exponential decay in the above has to be replaced by \\PnSn - PnSn\\2 < ^T - 3. Computation of Battle-Lemarié's wavelets Fix a positive integer m. Let A = B2m be the bi-infinite matrix whose symbol is Y^k€zB2m(k)e~Jko}. Clearly, A is a banded Toeplitz matrix. To see that A is a positive operator on I2, we show that A > ci for some c > 0 as follows: For any x e I2 , we have x'Ax=j-[ Fix)ico)FiA)i(ú)Fix)ioj) 2n J-n = FiA)i0^f dco \Fix)iœ)\2dœ >minFiA)i(ú)\\x\\22. w With c = minwFiA)iœ) > 0, we have A > ci. Similarly, we can show that License or copyright restrictions may apply to redistribution; see https://www.ams.org/journal-terms-of-use ON THE COMPUTATIONOF BATTLE-LEMARIÉ'SWAVELETS 697 11-4 - I\\2 < 1. Indeed, ||(4 - 7)x||2 = ¿ |* \FiA - I)ico)\2\F{x)iœ)\2 dco = ^Jjl-FiA)ico)\2\Fix)ioj)\2dco <max|l -FiA)ico)\2\\x\\22< ( 1 - minF(y4)(w) Nil Thus, we have ||M-/)x||2<(l-mjnF(^)(ta))||x|b and hence, \\A-1\\2 < 1. Thus, B2TOsatisfies all the conditions of Theorem 1. By (1), we have ., . 1 /sin&>/2\m y/EkezB2m(k)e-Jk- V co/2 J Thus, <PÍt)= ¿ZkakBmit-k) with x = iak)keZ satisfying C2mx = Ô and C22m= B2m . Using our Theorem 1, we conclude that our numerical method is valid to compute the ak\. By (2), the transfer function is H ico) = v -cosm(w/2). Note that when m is even, then cosm(w/2) = (1 - (ejo>+ e-J(°)/2)ml2, which is a finite series. However, when m is odd, cosm(w/2) is no longer a finite series. In order to compute i/(<y), let Sm be the Toeplitz matrix whose symbol is cos2m(<y/2) = (1 - iejw + e~Jw)/2)m . Let Z be a zero insertion operator on I2 defined by 7 7t \ i ^ Zx = Z(x/)/ez = (z/)/€Z :u Í x'/2 with z¡, = i ' i 0 if'is even, if ms odd. Thus, Hi(o) = Yikezßke-ikw with x = ißk)keZ satisfying x = w*y *z, where * denotes the convolution operator of two vectors in I2 and y = CmS, z = ZC~xS, w = T¿ with C2, = B2m , T2, = Sw . Using our Theorems 1 and 2, we know that our numerical method gives a good approximation to y and z. For m even, our numerical method produces an xn which converges pointwise to x exponentially. When m is odd, the remark after Theorem 2 has to be applied, and the Wat produced by this procedure does no longer converge to w exponentially. License or copyright restrictions may apply to redistribution; see https://www.ams.org/journal-terms-of-use 698 MING-JUN LAI By (3), the wavelet y/ associated with tp is given by y/{2co) = e-JO}Hico + n)$ito). Once {ak}k€Z and {ßk}kez are computed, {yk}kez can be obtained by convolution. We have implemented this method to compute Battle-Lemarié's wavelets in MATLAB. The graphs of Battle-Lemarié's wavelets are shown in the following figures. 2 ~i i i i i r T r 1.5 1 0.5 0 .0.5 -1 0 2 4 6 8 -10 -8 -6 -4 -2 0 Battle-Lemarié's scaling function and wavelet of degree 1 1.5 J-1—j_i_ 2 4 10 "i-1-1-1-r 1.0 0.5 -0.5 \-15 -10 0 scaling function and wavelet of degree 3 Battle-Lemarié's 1.2 i i i 1.5 r 1.0 10 15 ~\-1-1-r 1.0 0.8 0.6 0.5 0.4 0.2 0 -0.5 -0.2 - -20 -15 -10 -5 0 5 10 Battle-Lemarié's 15 20 -20 -15 -10 -5 0 scaling function and wavelet of degree 5 License or copyright restrictions may apply to redistribution; see https://www.ams.org/journal-terms-of-use 5 10 15 20 ON THE COMPUTATIONOF BATTLE-LEMARIÉ'S WAVELETS -25 -20 -15 -10 -5 0 5 10 15 20 25 Battle-Lemarié's -25 -20 -15 -10 -5 0 scaling function and wavelet of degree 7 699 5 10 15 20 25 1.5 1.0 0.5 -0.5 -25 -20 -15 -10 -5 0 5 10 15 20 25 Battle-Lemarié's -30 -20 -10 0 scaling function and wavelet of degree 9 10 20 30 Bibliography 1. G. Battle, A block spin construction of ondelettes, I: Lemarié functions, Comm. Math. Phys. 110(1987), 601-615. 2. I. Daubechies, Ten lectures on wavelets,SIAM, Philadelphia, PA, 1992. 3. P.-G. Lemarié, Ondelettes à localisation exponentielles, J. Math. Pures Appl. (9) 67 (1988), 227-236. 4. S. Mallat, Multiresolution approximations and wavelet orthonormal bases of L2(R), Trans. Amer. Math. Soc. 315 (1989), 69-87. 5. _, A theory for multiresolution signal decomposition: the wavelet representation, IEEE Trans. Pattern Anal, and Machine Intelligence 11 (1989), 674-693. Department of Mathematics, E-mail address : mj laiQw iener. University of Georgia, Athens, Georgia math. uga. edu License or copyright restrictions may apply to redistribution; see https://www.ams.org/journal-terms-of-use 30602