Hajji MA^{*}
Department of Mathematical Sciences, United Arab Emirates University, United Arab Emirates
Received Date: January 12, 2016; Accepted Date: February 01, 2016; Published Date: February 05, 2016
Citation: Hajji MA (2016) On the Exact Values of Daubechies Wavelets. J Phys Math 7:157. doi:10.4172/2090-0902.1000157
Copyright: © 2016 Hajji MA. This is an open-access article distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original author and source are credited.
Visit for more related articles at Journal of Physical Mathematics
In this paper we propose an algorithm for the calculation of the exact values of compactly-supported Daubechies wavelet functions. The algorithm is iterative, performing a single convolution operation at each step. It requires solving, at the first step only, a linear system of a relatively small size. The novelty of the algorithm is that once the values at dyadic points at a certain level j are calculated they do not need to be updated at the next step. We find that this algorithm is superior to the well-known cascade algorithm proposed by Ingrid Daubechies. This algorithm can serve well in wavelet based methods for the numerical solutions of differential equations. The algorithm is tested on Daubechies scaling functions as well as Daubechies coiflets. Comparison with the values obtained using the cascade algorithm is made. We found that the cascade algorithm results converge to ours
Daubechies wavelets; Daubechies coiflets; Multiresolution analysis; The cascade algorithm; Convolution
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].
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.
Wavelet basis is a doubly-index family of L^{2}(R) functions, ψ_{j,k}, j,k∈Z, defined by
(1)
where ψ(x) is the mother wavelet defined in terms of a mother scaling function φ(x) via the refinement equation:
(2)
The mother scaling function, φ(x), is itself defined recursively by the refinement equation:
(3)
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
• and .
•
•
• is an orthonormal basis of V_{0}, where
A consequence of the above MRA sructure are the following:
• The set of functions is an orthonormal basis for V_{j}, where .
• Associated with each V_{j} there is a wavelet space W_{j}, its orthogonal complement in V_{j+1}, i.e., .
• An orthonormal basis for W_{j} is the set .
•.
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 f_{j}∈ V_{j} through
Similarly, the orthogonal projection of f onto W_{j} gives
By the structure of the MRA , we have
.
The coefficients at a lower scale are computed from via Mallat’s algorithm [17]:
and conversely are reconstructed via
(6)
For compactly-supported wavelets, the sums in (2) and (3) are finite:
with 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 is an iterative scheme proposed by Daubechies explanied in [1-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
(9)
(10)
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,
(11)
Moreover if f is uniformly continuous, then the convergence in (11) is uniform as well, and If f is lder continuous with exponent a, i.e.,, then the convergence in (11) is exponentially fast in j, i.e.,
(12)
A consequence of the above proposition is the following.
Lemma 2 For any dyadic rational
(13)
Moreover, if j is sufficiently large, say we have the estimate
(14)
Where C and j_{0} depend on m and k.
Proof By Proposition 1, we have
which proves (13). Estimate (14) follows from (12).
Lemma 2 suggests that at any j -level, φ(2^{j} x) can be approximated by
(15)
For j,k∈Z , define the coefficients . Then by (15), we have
(16)
The coefficients can be reconstructed recursively using Mallat’s algorithm (6) starting at the scale j=0. At scale j=0, by (9), we have
(17)
Since for all It follows from (6) that are given by
(18)
The cascade algorithm summarizes as follows:
• Start with the sequence which can be viewed as a first approximation of at the integers.
• For recursively via (18). The sequence gives the approximation of φ at the j -level dyadics .
We note that if the length of the filter h is L then the length of the sequence is. At every step j , the algorithm computes for the first time approximations of φ at the odd dyadics and refines the approximations at the even dyadics which were obtained at the previous j-1 step. To be more precise at step j, the algorithm computes a total of values (the length c_{k}^{j} ) of which values comprise the new approximations (at the odd dyadics) and the rest values constitute refinements of the old approximations obtained at the previous j-1 step.
Before closing this section, we would like to mention that at any scale j the cascade algorithm does not cover all j -level dyadic points in the interval [0, L −1] ; it only covers dyadics for Since there are j-level dyadics in [0, L −1] , the (L − 3) dyadic points for k = 0,1,...., L − 4 , are not covered; some of them, but not all, will be covered at next step j +1 .
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 x = 2^{− j}m, j,m∈Z . 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 n =1,2,.., L − 2 (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, L −1] , we have, by continuity, φ(x) =0 for x ≤ 0 and x ≥ L −1 . Let
(19)
be the column vector containing the values of φ at the integers. Then, according to the refinement equation,
(20)
Φ^{(0)} satisfies the linear system
Φ^{(0)}=AΦ^{(0) } (21)
where A is an (L-2)×(L-2) matrix with entries , given by
(22)
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.
Proposition 3 Suppose that there exists a continuous solution φ to (20). Then λ=1 is an eigenvalue of the matrix A in (22) and Φ^{(0)} is an associated eigenvector.
Proof Let Φ^{(0)} be as in (19). Then by (20) Φ^{(0)} satisfies (21). If Φ^{(0)}=0 (the zero vector), then φ(n)=0 for . This implies, using (20), that . Then by continuity of φ, this would imply that φ=0 which is a contradiction. It follows that Φ^{(0)}≠0 and consequently Φ^{(0)} is an eigenvector of A associated to the eigenvalue λ=1.
Since φ satisfies for any x, . Then the vector Φ^{(0)} is equal to the normalized eigenvector of A corresponding to λ=1 . Precisely, Φ^{(0)} is the unique solution to the (L-1)×(L-2) nonhomogeneous system
Bx=b (23)
where
and 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
(24)
Note that we do need to calculate φ at the even half integers as they were computed in Φ^{(0)}. Let Φ^{(1)} be the column vector containing the values of φ at the odd 1-level dyadics Then (24) can be viewed as a convolution operation followed by downsampling by 2,
(25)
where “conv” denotes convolution and ↓ denotes downsampling by 2.
Let , be the vector of length containing the values of φ at the odd j-level dyadics, i.e.,
(26)
A careful examination of the refinement equation (7) gives
(27)
where is an “upsampled” version of the vector , obtained by inserting zeros between every two successive entries in . Explicitly,
(28)
The length of the vector is equal to Note that the length of ^{(j)} given by the convolution formula (27) is the sum of the lengths of less one, i.e.,, and it agrees with the length of Φ^{(j)} given by formula (26).
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 and downsample by 2 to get Φ^{(1)}. This gives the values of φ at the odd 1-level dyadics .
3. For j ≥ 2, compute the values of φ at the odd j-level dyadics by convolving of the vector with the vector Φ^{(j-1)} where, again, Φ^{(j-1)} is the vector containing only the values of φ at the odd (j-1)-evel dyadics.
The values of the wavelet function ψ(x) at the any j-level dyadics can be computed using the relation
(29)
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 φ, , the sought eigenvector Φ^{(0)} can be obtained as the solution of (23).
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 crucial one because it yields the exact values at the integers from which everything else is derived.
• 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.
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.
Values of Daubechies’ scaling function db4 | |||||
---|---|---|---|---|---|
* x | The cascade algorithm | Our algorithm | |||
j = 5 | j =10 | j =15 | j = 20 | j = 0 | |
1.00747364 | 1.00717932 | 1.00717027 | 1.00716998 | 1.00716997 | |
-0.03432410 | -0.03385159 | -0.03383741 | -0.03383696 | -0.03383695 | |
0.03983843 | 0.03961669 | 0.03961065 | 0.03961046 | 0.03961046 | |
-0.01180456 | -0.01176501 | -0.01176437 | -0.01176435 | -0.01176435 | |
-0.00120268 | -0.00119824 | -0.00119796 | -0.00119795 | -0.00119795 | |
0.00001926 | 0.00001883 | 0.00001882 | 0.00001882 | 0.00001882 |
Table 1: Values of db4 at the integers.
Values of Daubechies’ scaling function db6 | |||||
---|---|---|---|---|---|
*x | The cascade algorithm | Our algorithm | |||
j = 5 | j =10 | j =15 | j = 20 | j = 5 | |
0.0004331 | 0.0002749 | 0.0002707 | 0.0002705 | 0.0002705 | |
0.1696324 | 0.1623220 | 0.1620967 | 0.1620896 | 0.1620894 | |
0.8333298 | 0.8208984 | 0.8205067 | 0.8204944 | 0.8204941 | |
0.9002129 | 0.9135170 | 0.9139223 | 0.9139350 | 0.9139354 | |
-0.1721720 | -0.1591060 | -0.1586933 | -0.1586804 | -0.1586800 | |
0.1241885 | 0.1214656 | 0.1213838 | 0.1213812 | 0.1213811 | |
-0.0268479 | -0.0277469 | -0.0277748 | -0.0277757 | -0.0277757 | |
0.0014295 | 0.0014226 | 0.0014213 | 0.0014212 | 0.0014212 | |
-0.0022161 | -0.0023020 | -0.0023046 | -0.0023047 | -0.0023047 | |
0.0000263 | 0.0000270 | 0.0000270 | 0.0000270 | 0.0000270 |
Table 2: Values of db6 at selected level dyadic.
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.
Make the best use of Scientific Research and information from our 700 + peer reviewed, Open Access Journals