6

THE DAUBECHIES WAVELETS

The wavelets that we have looked at so far—Haar, Shannon, and linear spline wavelets—all have major drawbacks. Haar wavelets have compact support, but are discontinuous. Shannon wavelets are very smooth, but extend throughout the whole real line, and at infinity they decay very slowly. Linear spline wavelets are continuous, but the orthogonal scaling function and associated wavelet, like the Shannon wavelets, have infinite support; they do, however, decay rapidly at infinity.

These wavelets, together with a few others having similar properties, were the only ones available before Ingrid Daubechies discovered the hierarchy of wavelets that are named after her. The simplest of these is just the Haar wavelet, which is the only discontinuous one. The other wavelets in the hierarchy are both compactly supported and continuous. Better still, by going up the hierarchy, they become increasingly smooth; that is, they can have a prescribed number of continuous derivatives. The wavelet’s smoothness can be chosen to satisfy conditions for a particular application. We now turn to Daubechies’ construction of the first wavelet past Haar.

6.1 DAUBECHIES’ CONSTRUCTION

Theorem 5.23 lists the three sufficient conditions on the polynomial P that ensures that the iteration given in Section 5.3.4 produces a scaling function.

For a given polynomial P(z), let

In terms of the of the function p, the three conditions in the hypothesis of Theorem 5.23 can be stated as

(6.1)

(6.2)

(6.3)

In this section, we describe one polynomial, due to Daubechies, which satisfies (6.1)–(6.3).

From Example 5.22, the polynomial associated with the Haar scaling function is

This choice of p₀ satisfies Eqs. (6.1)-(6.3). However, the Haar scaling function is discontinuous. One way to generate a continuous scaling function is to take convolution powers. In fact, the convolution of the Haar scaling function with itself can be shown to equal the following continuous function (see exercise 5 in Chapter 2):

(see Figure 6.1). Now the Fourier transform of a convolution equals the product of the Fourier transforms. In particular, the Fourier transform of ϕ * ϕ * … ϕ (n times) is In view of Eq. (5.29) for the Fourier transform of the scaling function, we may be first led to try for some suitable power of n, as the function that generates a continuous scaling function. However, property (6.2) no longer holds (unless n = 1, the Haar cases.

Instead of simply raising p0 to the nth power, we raise both sides of the identity to the nth power. With n = 3, we obtain

(6.4)

or

Using the identities cos(u) = sin(u + π/2) and sin(u) = – cos(u + π/2) for the last two terms on the right, we have

If we let

then the previous equation becomes

and property (6.2) is satisfied. Property (6.3) is also satisfied since Also note that |p(0)| = 1. So all that is left is to identify p itself (we have only defined |p|). First, we rewrite the defining equation for |p| as

By taking square roots of this equation, we let

where a() is a complex-valued expression with |a()| = 1 which will be chosen later.

To identify the polynomial P (with ), we use the identities

to obtain

We choose in order to clear all positive and fractional powers in the exponent. Expanding out and collecting terms, we obtain

The equation is therefore satisfied by the polynomial

Since p satisfies (6.1) through (6.3), P satisfies the hypothesis of Theorem 5.23. Recall that Therefore, for Daubechies’ example, we have

(6.5)

The Daubechies scaling function, ϕ, is found by applying the iterative procedur described in Theorem 5.23. Figures 5.18 through 5.21 in Chapter 5 illustrate th first few iterations of this procedure. Figure 6.2 shows the (approximate) graph of ϕ that results from iterating this procedure many times.

What about the wavelet ψ that is associated with the Daubechies scaling function ϕ? As shown in Theorem 5.10, once the coefficients, p_k, have been identified and ϕ has been constructed, then the associated wavelet is given by the formula

Figure 6.3 shows the (approximate) graph of the associated wavelet function.

Unlike the Haar scaling and wavelet functions, the Daubechies scaling and wavelet functions are continuous. However, they are not differentiable.

6.2 CLASSIFICATION, MOMENTS, AND SMOOTHNESS

Other, smoother scaling functions and wavelets can be obtained by choosing a higher power than n = 3 in Eq. (6.4). Any odd power, n = 2N - 1, can be used there. In fact, Daubechies showed that for every N there will be 2N nonzero, real scaling coefficients p0, . . ., p_2n-1, resulting in a scaling function and wavelet that are supported on the interval 0 ≤ t ≤ 2N – 1. They are chosen so that the corresponding degree 2N – 1 polynomial has the fac-torization,

(6.6)

where the degree of is N – 1 and This guarantees that the associated wavelets will have precisely N “vanishing moments.” We will discuss what this means and why it is important later.

Apart from a reversal of the coefficients in P_N (cf. exercise 1), these coefficients are unique. In the two cases that we have dealt with, N = 1 (Haar) and N = 2 (Daubechies), the polynomials are and P2(z) = Both polynomials have the factorization in Eq. (6.6),

with and

The scaling function ϕN and wavelet ψ_N generated by P_N have Fourier transforms given in terms of infinite products. Because the coefficients of P_N are real, Thus from Eqs. (5.29) and Eq. (5.30), we have that

(6.7)

(6.8)

Note that because P_N(–1) = 0. If N > 1, we also have because In general, by exercise 2, we have that

(6.9)

This gives the following result.

Proposition 6.1 For the Daubechies wavelet ψ_N, we have

(6.10)

Proof. The proposition follows from Eq. (6.9) and the Fourier transform property

Just set f = ifrff, n = k, and k = 0.

In mechanics, integrals of the form called the moments of a mass distribution p. The term “moment” carries over to an integral of any function against x^k. We can thus rephrase the proposition by saying that ψNhas its first N moments vanish. This is usually shortened to saying that ψN has N vanishing moments; that they are the first N is understood.

Daubechies wavelets are classified according to the number of vanishing moments they have. The smoothness of the scaling function and wavelet increase with the number of vanishing moments. The N = 1 case is the same as the Haar case; both scaling function and wavelet are discontinuous. The N = 2 Daubechies scaling function and wavelet are continuous, but certainly do not have smooth derivatives. In the N = 3 case, both are continuously differentiable. When N is large, the number of continuous derivatives that ϕN an ψN have is roughly N/5. So, to get 10 derivatives, we need to take N ≈ 50. In the following table, we list approximate scaling coefficients for the Daubechies wavelets, with N going from 1 to 4. Of course, the scaling coefficients given for N = 2 are just the decimal approximations for those found in Eq. (6.5).

Why is it useful to have vanishing moments? The short answer is that vanishing moments are a key factor in many wavelet applications: compression, noise removal, and singularity detection, for example.

To understand this better, let’s look closely at the N = 2 case. According to Eq. (6.10) with N = 2, the first two moments (k = 0,1) vanish. The k = 2 moment is From our earlier discussion, we see that for the third (k = 2) moment, we have

with

(6.11)

We want to use these moments to approximate the wavelet coefficients for smooth signals, and we want to show that these coefficients will be small when the level j is high. If f is a smooth, twice continuously differentiable signal,

then its j, k-wavelet coefficient is

If j is large enough, the interval over which we are integrating will be small, and we may then replace by its quadratic Taylor polynomial in Doing this, we get the following approximation for b^k_j :

(6.12)

The integral on the right can be reduced to doing the integrals for the first three moments of ϕ. Since the first two moments vanish, and the third is given in (6.11), we can evaluate the integral. (See exercise 3.) Our final result is that the j, k-wavelet coefficient is approximately

(6.13)

Singularity Detection. As an application of the formula in (6.13), we find a point where an otherwise smooth function has a discontinuity in its derivative. This is called singularity detection, and this process can be used, among other things, to detect cracks in material.

As an example, we will determine where the piecewise linear function shown in Figure 6.4 changes slope. Keep in mind that the formula (6.13) is exact in any region where a signal f is constant, linear, or quadratic in x, because the quadratic Taylor polynomial for f is f in these cases. Since where f is linear, formula (6.13) implies that the only nonzero wavelet coefficients will come from a small region near the corner point where the slope changes. The signal itself has the form

After sampling it at 256 equally spaced points, we decompose it using the N = 2 Daubechies wavelet. Thus, our starting level was 7=8; the level down from this is j = 7. The only appreciable wavelet coefficient was The rest were on the order of 10^–14. The index k = 97 corresponds to a wavelet supported on the x interval [2.79, 2.83]; consequently, the singularity is located in this interval.

6.3 COMPUTATIONAL ISSUES

Consider the following problem. We want to use N = 2-Daubechies’ wavelets to decompose a signal that we have n samples for, s₀ …, s_n–1. These samples will be regarded as part of our top-level approximation coefficient sequence, aⁱ. We will filter these numbers using the high- and low-pass filters for the Daubechies wavelets. These are H and l, and they have impulse responses given by the sequences h and I, where

We convolve h and h with a⁷ and down sample. Thus, for any a^j—not just our signal with eight samples—we have

(6.14)

(6.15)

(We use k + 1 in (6.15) rather than k in order to deal with the same a^j’s on the right side in both formulas.)

To compute each level j – 1 coefficient, we need four consecutive samples, always starting at an even number. For example, if we have n = 8 samples, then to compute the a₂² coefficient, we need s₄ through s₇. To compute a₃², we need s₆ through s₉. But we do not have s₉! It’s no surprise that at some point we will run out of coefficients to filter. The surprise is that we are only able to compute k = 0, 1, 2. This means from eight samples, we only get three decomposition coefficients, not the four that we would expect. This is the overspill problem. The filters H and L spill over to samples we don’t have.

What to do? First of all, the overspill problem comes about because we don’t know what the signal was before or after our set of samples. Thus we need to extend the signal in some fashion beyond our initial set of samples. Here are a few of the standard choices. In the accompanying illustrations, the filled circles represent the original signal, and open circles are points in the signal’s extension.

Zero-Padding. Extend the signal by adding zeros at the two ends. This is appropriate if the signal is very long, and the ends of the signal do not matter, or if the signal really is abruptly started and stopped. This amounts to setting s_k = 0 if k π 0 or k > n-1.

Periodic Extension. Another approach is to re-use the sampled data by making the signal a periodic one, so that s_k+n = s_k. For example, with eight samples, s₀, … , s₇, we would let s₈ = s₀, s₉ = s₁, and so on.

Smooth Padding. Extend the signal by linearly extrapolating the data near the two ends. This is appropriate if signal isn’t too noisy, at least near the ends.

Symmetric Extensions. The signal is evenly extended at the endpoints by reflection. There are two ways to do this. One can reflect about a line through an endpoint, or about a line midway between an endpoint and the next point on the grid. The first type is illustrated at the left endpoint, and the second type is shown at the right endpoint.

6.4 THE SCALING FUNCTION AT DYADIC POINTS

Although the values of the scaling function ϕ do not enter into the decomposition and reconstruction algorithms, it is useful to be able to compute approximate values of ϕ in order to verify some of its properties (like continuity). An iterative method for computing ϕ is given in Theorem 5.23. However, this algorithm is somewhat cumbersome from a computational point of view. A more computationally efficient method for computing the scaling function ϕ is to use the scaling equation to compute the value of ϕ at all dyadic values, x = 1/2ⁿ, where I and n are integers. We illustrate this process in the following steps. To simplify the exposition, we shall concentrate on the Daubechies N = 2 case with four nonzero p coefficients as constructed in the previous section, but the procedure easily generalizes.

Step 1. Compute = at All the Integer Values. Let ϕ_l = ϕ(l) for l ∈ Z. The Daubechies (N = 2) scaling function is nonzero only on the interval 0 < x <3 with ϕ₀ = ϕ₃ = 0 (see Figure 6.2). In particular, ϕ₁ and ϕ2 are the only nonzero unknown values of ϕ at integer points. In order to arrange the normalization ϕ = 1, we require or in our specific case:

(6.16)

Now we use the scaling equation For x = l, this equation becomes

At x = 2, the scaling equation becomes

These two equations can be written in matrix form as

Here, the p values are known [see Eq. (6.5)] and we are solving for the unknowns ϕ₁ and ϕ₂. In order for this matrix equation to have a nonzero solution, the matrix must have an eigenvalue equal to 1. Then (ϕ₁,ϕ₂) would be the corresponding eigenvector with ϕ₁ + ϕ₂ = 1. To find this eigenvector, we rewrite the matrix equation as

A nonzero solution to this matrix equation exists if the first row is a multiple of the second. From Theorem 5.9, ∑ _Podd = ∑ _Peven = 1. Therefore the first row is the negative of the second. The equation corresponding to the first row is

We restate the normalization equation (6.16):

These two equations can be solved simultaneously. In the case of Daubechies, the solution is

The ϕ values at all other integer points are zero. Therefore, all the ϕ_t = ϕ(l), l ∈ Z, have been determined.

Step 2. The Values of ϕ at the Half-Integers. To compute ϕ(l/2), we use the scaling equation

(6.17)

at x = l/2 to obtain

(6.18)

Since ϕ(l – k) is known from Step 1, ϕ(l/2) can be computed. We need only compute ϕ(l/2) for l = 1, 2, 3,4, and 5 since ϕ(x) = 0 for x π 0 and x > 3. When l = 2 or 4, l/2 is an integer and the values of ϕ at these points are known from Step 1. Thus we need only compute (2) for l = 1, 3, and 5. Equation (6.18) for / = 1, 3, 5 becomes

The condition ∑_lϕ (l) = 1 (from Step 1) implies ∑_lϕ (l/2) = 2. Indeed, Eq. (6.17) gives

By the change of index, the inner sum is This sum equals 1 from Step 1. By Theorem 5.9, ∑_kP_k = 2. Therefore

as claimed.

Step 3. Iteration. The values of ϕ at the quarter integers, l/4, can be computed from the values of ϕ at the half-integers by letting x = l/4 in the scaling equation In general, once ϕ has been computed at the values x = l/2^n-1, then we can compute the value of ϕ at the values x = l/2ⁿ by inserting x = l/2ⁿ into the scaling equation (6.17):

The right side involves values of ϕat which have been computed from the previous step.

We claim

(6.19)

We have already shown this equation for n = 0 (Step 1) and n = 1 (Step 2). Suppose by induction that we assume this equation is true for n – 1. We now show it to be true for n. We have

By the induction hypothesis, the inner sum is 2^n-1 [Eq. (6.19) with n replaced by n – 1]. Therefore

Since ∑_k Pk = 2 (Theorem 5.9), the right side equals 2ⁿ", as desired.

As n gets larger, the set of dyadic points {l/2ⁿ, l ∈ Z} get denser. Since any real number is a limit of dyadic points and since the Daubechies scaling function is continuous, the value of ϕ at any value x can be obtained as a limit of ϕ-values at Dyadic points.

The scaling function, ϕ that is constructed in this manner, satisfies the normalization condition ∫ ϕ = 1. To see this, we regard ∫ ϕ dx as a limit as n → ∞ of a Riemann sum over a partition given by {x_t = l/2ⁿ; l = … – 1,0,1, … } The width of this partition is Δx = 1/2ⁿ. So, we have

Since the right side is 1, as claimed.

EXERCISES

1. Show that satisfies the N = 2 Daubechies scaling relation, but with the coefficients reversed; that is,

What is the corresponding equation for Sketch both and

2. Show that Eq. (6.9) holds.

3. Use the first three moments of Ψ = Ψ_N-1=2 to evaluate integral in Eq. (6.12) and thus obtain the approximation to b^J_k given in Eq. (6.13).

4. Repeat exercises 9 and 10 in Chapter 4 using Daubechies N = 2 wavelets.

5. Repeat exercise 4 for the signal

where

Compare your results using Daubechies wavelets with those of exercise 9 in Chapter 3.

6. Let f (t) be defined via

Sample f at the dyadic points k × 2^-8, k = -256 … 512. Thus, j =7 is the top level. Using the Daubechies N = 2 wavelet, implement a one-level decomposition. Plot the magnitudes of the level 7 wavelet detail coefficients. Which wavelet has the largest coefficient? What t corresponds to this wavelet? Explain.

7. Let g be defined by

Sample g at the dyadic points k × 2^-8, k = -256... 512. (j = 7 is the top level.) Using the Daubechies N = 2 wavelets, implement a one-level decomposition. Plot the magnitudes of the level 7 wavelet coefficients for each, and determine which wavelet coefficient is greatest. Repeat with the N = 3 wavelet. Give an explanation of why the greatest coefficients occur where they do.

8.(Polynomial Suppression) Frequently, signals have a “bias” that is a polynomial part (usually linear or quadratic) added to a bounded rapidly oscillating part–for example, Suppose that we have taken 1024 samples of f on the interval [-1,1]. Come up with a strategy for separating the two signals via a Daubechies wavelet analysis. Use MATLAB to carry it out. (Hint: What is the smallest value of N needed to reproduce quadratics exactly? Also smooth padding works best for these examples.)

9.(Noise Removal) Wavelets can be used to remove noise from a signal. Let f(t) = sin(8πt) cos(3πt) + n(t), with n(t) being the noise. Numerically, we can model n(t) by using a random number generator, such as MATLAB’s rand. Take 1500 samples on [-2, 3] of f and do an analysis with N = 2,3, and 6 Daubechies wavelets. Experiment with different levels. Which wavelet does the best job?