On the Exact Values of Daubechies Wavelets

The widespread interest in wavelets and their applications started in the 1980s after the breakthrough made by Daubechies [1,2] in constructing the first orthogonal compactly-supported wavelets with arbitrary regularity. Since then many researchers from different fields of science and engineering jumped into the world of wavelets with different intentions. Some were interested in ways to apply wavelets in their fields and others were interested in developing new theories and generalizations.In the engineering side, wavelets have found great success in signal processing such as in the analysis of sound patterns and image processing [3-6]. Wavelets have also found great success in the design and efficient implementation of numerical algorithms for the solution of differential equations [7-16].


Introduction
The widespread interest in wavelets and their applications started in the 1980s after the breakthrough made by Daubechies [1,2] in constructing the first orthogonal compactly-supported wavelets with arbitrary regularity. Since then many researchers from different fields of science and engineering jumped into the world of wavelets with different intentions. Some were interested in ways to apply wavelets in their fields and others were interested in developing new theories and generalizations.In the engineering side, wavelets have found great success in signal processing such as in the analysis of sound patterns and image processing [3][4][5][6]. Wavelets have also found great success in the design and efficient implementation of numerical algorithms for the solution of differential equations [7][8][9][10][11][12][13][14][15][16].
Wavelet collocation based methods for solving different equations require knowledge of the values of the wavelet basis elements at collocation points. However, many of the available wavelets, such as the well-known Daubechies compactly-supported wavelets, do not have explicit formulas. Instead, they are defined recursively by refinement equations. In many applications having accurate values of the wavelet bases functions is very important in obtaining accurate solution to the problem. In [2] Daubechies described an algorithm known as the cascade algorithm for computing approximate values of the compactly-supported scaling and wavelet functions with arbitrary high precision. This algorithm works as a refinement scheme. At each step approximately twice as many values are computed, values at odd dyadic points 2 -j (2k+1) are computed for the first time and values at even dyadic points 2 -j (2k) are refined from the previous step. In the long run, i.e., as j→∞, the cascade algorithm produces the exact values. In this work, we propose an algorithm to calculate the exact values of Daubechies scaling and wavelet functions. The proposed algorithm avoids the refinement step in the cascade algorithm. Moreover, our algorithm, at each step computes only values at odd dyadic points 2 -j (2k+1). The values at even dyadic points 2 -j (2k) have already been computed at the previous step and no need for refinement because the values are exact.
The paper is organized as follows. In section 2, we briefly review compactly-supported scaling and wavelets functions and some related properties. In section 3, we outline the cascade algorithm in [2]. In section 4, we describe the proposed algorithm. In section 5, we apply the proposed algorithm to Daubechies scaling functions and compare our results to those obtained by the cascade algorithm. Finally, we conclude by some remarks in section 6.

Preliminaries
Wavelet basis is a doubly-index family of L 2 (R) functions, ψ j,k, j,k∈Z, defined by where ψ(x) is the mother wavelet defined in terms of a mother scaling function φ(x) via the refinement equation: The mother scaling function, φ(x), is itself defined recursively by the refinement equation: The coefficients h k and g k are known, in the language of signal processing, as the low-and high-pass filter coefficients, respectively. For orthogonal wavelets, they are related by The scaling function φ generates an orthogonal multiresolution analysis (MRA) [17] which is an increasing sequence of subspaces V j of L 2 (R) (approximation spaces) with the following properties , A consequence of the above MRA sructure are the following: • The set of functions { j,k (x), k∈Z} is an orthonormal basis for V j , where φ j,k (x)=2 j/2 φ(2 j x-k).
• Associated with each V j there is a wavelet space W j , its orthogonal complement in V j+1 , i.e., V j ⊕ W j = W j+1 .
• An orthonormal basis for W j is the set {ψ j,k (x), k∈Z}.
Let P j and Q j be the orthogonal projections onto V j and W j , respectively. Any function f ∈ L 2 (R) can be approximated by a function Similarly, the orthogonal projection of f onto W j gives , , , By the structure of the MRA ( and conversely j k s are reconstructed via For compactly-supported wavelets, the sums in (2) and (3) are finite: so that both and ψ are supported in [0, L-1] .The integer L=2M where M is the number of vanishing moments of ψ: [2]

The Cascade Algorithm
The cascade algorithm is an iterative scheme proposed by Daubechies explanied in [1][2][3] to calculate approximate values of the scaling and wavelet functions, φ and ψ, at rational dyadic points x=2 -m k . In order to compare this algorithm to our proposed algorithm, it is worthwhile to describe it.
The cascade algorithm is based on the key fact that the scaling function φ is the unique function satisfying 0, It is also based on the fact that 2 j φ(2 j x) is an approximate δ-function as j→∞ in the sense of the following proposition [2].

Proposition 1 If f is a continuous function. Then for any x∈R,
Moreover if f is uniformly continuous, then the convergence in (11) is uniform as well, and If f is H o  lder continuous with exponent a, i.e., | ( ) ( ) | | | , , , then the convergence in (11) is exponentially fast in j, i.e., A consequence of the above proposition is the following.

Lemma 2
For any dyadic rational x=2 -m k, Where C and j 0 depend on m and k.
Lemma 2 suggests that at any j -level, φ(2 j x) can be approximated by with the error (14)).
For j,k∈Z , define the coefficients . Then by (15), we have The coefficients j k c can be reconstructed recursively using Mallat's algorithm (6) starting at the scale j=0. At scale j=0, by (9), we have . It follows from (6) that j n c are given by 1 2 = .
The cascade algorithm summarizes as follows: • Start with the sequence 0 ,0 = n n c δ which can be viewed as a first approximation of at the integers.

The Proposed Algorithm
The algorithm we propose in this paper covers all dyadics at any given scale j. It also yields the exact values of the scaling and wavelet functions φ and ψ at dyadic points = 2 , , Since ψ is given in terms of φ, it is suffices to describe the algorithm for . The algorithm is based on finding the exact values of φ at the integers = 1, 2, , 2 n L −  (the 0 -level dyadics). These are found by solving an eigenvalue problem.
Once the values of φ at the integers are obtained, at each subsequent step, the algorithm performs a single convolution operation to calculate the values of φ at only odd-dyadics. The algorithm is explained in the following.
Since φ is supported in [0,1] L − , we have, by continuity, φ(x) =0 for be the column vector containing the values of φ at the integers. Then, according to the refinement equation, Φ (0) satisfies the linear system where A is an (L-2)×(L-2) matrix with entries a ij , 1 , Equation (21) suggests that Φ (0) is an eigenvector of A corresponding to the eigenvalue =1 . The existence of the eigenvalue λ=1 is justified by the following proposition.
Since φ satisfies and I (L-2) is the identity matrix of order (L-2) .
Next, once the values of φ at the integers are known, we apply the refinement equation to find the values of φ at the odd half integers can be viewed as a convolution operation followed by downsampling by 2, where "conv" denotes convolution and ↓ denotes downsampling by 2.
Let Φ (j) 1 j ≥ , be the vector of length 2 j-1 (L-1) containing the values of φ at the odd j-level dyadics, i.e., A careful examination of the refinement equation (7) gives where j h  is an "upsampled" version of the vector 2 h , obtained by inserting (2 j-2 -1) zeros between every two successive entries in 2 h. Explicitly, The length of the vector j h  is equal to 2 j-2 (L-1)+1. Note that the length of (j) given by the convolution formula (27) is the sum of the lengths of j h  and Φ (j-1) less one, i.e., Our proposed algorithm summarizes in the following steps: 1. Compute the values of φ at the integers(x=1, 2,…,L-2)given by the normalized eigenvector of A corresponding to the eigenvalue 1 and collect them in a vector Φ (0) .
2. Convolve Φ (0) with the vector 2 h and downsample by 2 to get Φ (1) . This gives the values of φ at the odd 1-level dyadics The values of the wavelet function ψ(x) at the any j-level dyadics can be computed using the relation All of the steps in the proposed algorithm are computationally trivial except perhaps for the first one, where an eigenvector of A needs to be found. The matrix A being sparse, however, this is not very difficult. Using the normalization of φ, As a comparison between the two algorithms we note the following: • Our proposed algorithm is similar to the well-known cascade algorithm in that the computations of both algorithms involve convolutions, except for the first step in our proposed algorithm, where an relatively small size linear system has to be solved.
• The first step in our algorithm is the most expensive but it is the • The cascade algorithm begins by making an initial guess for φ(x) (a Dirac delta function), whereas our algorithm begins by computing φ(x) at the integers (0-level dyadics).
• A clear advantage of our algorithm is that (i) it provides the exact values and (ii) once the initial step has been performed, at every subsequent step, only the values of φ(x) at the odd dyadics need to be computed, the even dyadics at step j being the odd dyadics at step j-1. In contrast, the cascade algorithm requires refinement of φ(x) at the even dyadics.

Numerical Results
We have implemented the proposed algorithm in Matlab and tested it to produce the values of Daubechies' scaling functions as well as Daubechies' coiflets. Plots of Daubechies scaling functions db4 and db6 (M=4 and M=6 ) as well as Daubechies' coiflets coif2 and coif4 (M=2 and M=3), obtained by our algorithm, are displayed in Figures  1 and 2.
We compared our numerical values with the ones obtained using the cascade algorithm which is implemented in Matlab under the function "wavefun". Samples of the numerical values obtained are displayed in Tables 1 and 2. Table 1 displays the values of db4 at the integers obtained by the cascade algorithm for j=5,10,15, 20 and the ones obtained by our algorithm using only j=0, i.e., the solution of the system (23). Table 2 displays selected values of db6 at the j=5 level dyadics obtained by the cascade algorithm for j=5,10,15, 20 and the ones obtained by our algorithm using j=5. The results clearly show that the cascade algorithm results converge to ours as j tends to infinity.We remark that one has to iterate the cascade algorithm for larger value of j to obtain accurate results at a lower k level dyadics. For instance, from Table 2, we needed to iterate the cascade algorithm until j=20 to obtain as closer results to the ones obtained by our algorithm with only j=5.

Conclusion
In this paper, we have presented an efficient algorithm for the computation of the exact values of refinable functions, in particular Daubechies' scaling and wavelet functions. Our motivation for this work stems from the need for more accurate point values of the widely used Daubechies wavelets. This certainly will be useful in the numerical solutions of differential equations where wavelets are being used. Our proposed algorithm produces the exact values whereas the wellknown cascade algorithm produces approximate values. What is good about the proposed algorithm is that, at each step, it computes values only at odd dyadics. The cascade algorithm, however, at each step calculates new values and refines old ones. The only expensive step in our algorithm is the first one where we need to solve a relatively small size linear system. This first step is the crucial one in that it provides us with the exact values (to machine precision) at the integers from which the rest is derived.We believe that having exactly values of wavelet functions will give better results in wavelet based numerical schemes for the solution of differential equations.