ABOUT WAVELETS ON THE INTERVAL*

JEAN CHENG^{1} AND WEI-CHANG SHANN^{2}

Abstract. Wavelets on the closed interval are constructed step by step from the Daubechies’ wavelets with compact supports. We follow mostly the guidelines by Andersson, Hall, Jawerth and Peters, with modifications on the construction of “edge” wavelets. Then we construct the filter operators discussed by Cohen, Daubechies and Vial. A filtering process is suggested for two-dimensional discrete transforms based on wavelets on the interval.

Key words. wavelets on the interval, filter

1. Introduction. Applying the general theory of multiresolution approximation for L^{2}(ℝ) by Mallat [8] and Meyer, Daubechies [5] constructed a family of wavelets with compact supports. We index this family of wavelets by a parameter p. The mother wavelet ψ(x) has a support of length 2p − 1 and it is L^{2}-orthogonal to polynomials of degree ≤ p − 1. For wavelet theory and its applications in some branches of mathematics, we refer to Meyer [9].

The construction of various families of wavelets, their mathematical properties, their relations and applications to other disciplines are found in the “ten lectures” by Daubechies [6]. A concise overview of the compactly supported wavelets are found in the article by Strang [12]. Over there, the construction of this family of wavelets all start from a suitable sequence of scaling coefficients c_{k}, where c_{k} ≠ 0 for 0 ≤ k ≤ 2p − 1 and c_{k} = 0 otherwise. These c_{k} shall satisfy certain conditions such that they are uniqely determined and there is a scaling function ϕ(x) ∈ L^{2}(ℝ) such that

(1.1)  ϕ(x) = ∑_{k} c_{k}ϕ(2x − k)

with suppϕ(x) = [0, 2p − 1]. This ϕ(x) generates a multiresolution approximation V_{j} of L^{2}(ℝ). Let ϕ_{jk}(x) = \sqrt{2^{j}}ϕ(2^{j}x − k), then \{ϕ_{jk}(x) | k ∈ ℤ\} is an orthonormal basis of V_{j}. The mother wavelet is then defined by

ψ(x) = ∑_{k} g_{k}ϕ(2x − k).

Let ψ_{jk}(x) = \sqrt{2^{j}}ψ(2^{j}x − k), then \{ψ_{jk}(x) | (j, k) ∈ ℤ^{2}\} is an orthonormal basis of L^{2}(ℝ). Here g_{k} = (−1)^{k}c_{1−k}, c_{k} and g_{k} is a pair of perfect reconstruction quadrature mirror filters.

A construction of wavelets on the interval was first proposed by Meyer [10]. However, there are a few shortcomings in this construction. Two groups (see next paragraph) of researchers were carrying out alternative constructions with better properties and more stable computational procedures. They lead to very similar constructions, thus a joint announcement by both groups was written in [3].

In Section 2, we following the guideline by Andersson, Hall, Jawerth and Peters [1] to construct the scaling functions ϕ_{jk}(x) and the multiresolution approximation on L^{2}([0, 1]). In Section 3 we modify the approach in [1] to construct the wavelets ψ_{jk}(x) on L^{2}([0, 1]). This modification makes it more like the construction by Cohen, Daubechies and Vial [4]. In Section 4, again following [4], we describe the idea and construction of prefilters, then we propose a filtering process for twodimensional discrete wavelet transforms based on wavelets on the interval.

* This work was supported by the National Science Coucil in Taiwan; grant #84-2121-M-008-009.

1 E-mail address: zhenj@math.ncu.edu.tw.

2 Department of Mathematics, National Central University, Chung-Li, Taiwan, R.O.C. E-mail address:

shann@math.ncu.edu.tw

In the following, we sketch the construction step by step for the case of p = 2. At proper places we will indicate how would the parameters depend upon p. Note that all parameters can be computed from the scaling coefficients c_{k} and a few quadratures. Numerical methods are known for the computation of these quadratures in high precision [11].

2. Multiresolution on the interval. Let us assume the interval is [0, 1]. Near the left boundary x = 0, we define

X_{j}^{m}(x) = \sqrt{2^{j}}(2^{j}x)^{m};

near the right boundary x = 1,

Y_{j}^{m}(x) = X_{j}^{m}(x − 1) = \sqrt{2^{j}}(2^{j}(x − 1))^{m}.

Here m is an integer with 0 ≤ m ≤ (1 = p − 1).

We find the L^{2}-projection of X_{j}^{m}(x) in

V_{j}^{L} = span{ϕ_{jk}(x) | k = (−2 = 2 − 2p), −1, …, (0 = δ_{L} − 1)};

and Y_{j}^{m}(x) in

V_{j}^{R} = span{ϕ_{jk}(x) | k = (2^{j} − 3 = 2^{j} − 2p + 2 − δ_{R}), 2^{j} − 2, …, 2^{j} − 1}.

Here δ_{L} and δ_{R} are the parameters defined in [1] which indicate the number of “interior” basis functions that are included for the orthonormalization process near the boundary points. Following [4], we see that it makes more sense for δ_{L} ≥ 1 and δ_{R} ≥ 1. In this article we assume that δ_{L} = δ_{R} = 1 (same as in [4]). See remarks by the end of this section. We further assume that

j ≥ 3 = ⌈log_{2} 2p − 1⌉ + 1;

that is, V_{j}^{L} ∩ V_{j}^{R} = ∅.

Let x_{k}^{m} be the component of X_{j}^{m} in ϕ_{jk}(x):

x_{k}^{m} = (X_{j}^{m}, ϕ_{jk}) = ∫ (x + k)^{m}ϕ(x) dx = ∑_{i=0}^{m} \binom{m}{i}k^{i}M_{m−i}

Here M_{i} = ∫ x^{i}ϕ(x) dx, it is well known that these quadratures can be obtained as accurate as the scaling coefficients c_{k}. We see that x_{k}^{m} are independent of j. The components y_{k}^{m} are similarly defined and are independent of j.

Now we have the projected functions on [0, 1]:

X_{j}^{m} = ∑_{k=−2}^{0} x_{k}^{m}ϕ_{jk}|_{[0,1]}, Y_{j}^{m} = ∑_{k=2^{j}−3}^{2^{j}−1} y_{k}^{m}ϕ_{jk}|_{[0,1]}.

Obviously {X_{j}^{m}}_{m=0}^{1} and {Y_{j}^{m}}_{m=0}^{1} are linearly independent, and they are orthogonal with the “interior” scaling functions {ϕ_{jk}}_{k=1}^{2^{j}−4}, which are already orthonormal. See Figure 1, note that they are actually polynomials near the boundary points. So, we would only have to orthonormalize X_{j}^{m} and Y_{j}^{m} near the boundary points. The result will be an orthonormal set. Let V_{j}[0, 1] be the linear space spanned by this orthonormal set. Then

dim V_{j}[0, 1] = 2^{j}, j ≥ 3.

The orthonormalization process is done by taking the squared root inverse of the mass matrices of {X_{j}^{m}} and {Y_{j}^{m}}. Let X = (X_{mn}) where

X_{mn} = (X_{j}^{m}, X_{j}^{n}).

X is symmetric and positive definite, and it is actually independent of j:

\begin{aligned}X_{00} & = 1 + N_{1} + 2N_{2} + N_{3}, \\ X_{01} & = x_{−2}^{1}(N_{1} + N_{2}) + x_{−1}^{1}(N_{2} + N_{3}) + x_{0}^{1}, \\ X_{11} & = (x_{−2}^{1})^{2}N_{1} + (x_{−1}^{1})^{2}N_{3} + (x_{0}^{1})^{2} + 2x_{−2}^{1}x_{−1}^{1}N_{2} \end{aligned}

The accuracy of X is totally determined by the numerical quadratures

\begin{aligned}N_{1} & = ∫_{0}^{∞} ϕ(x + 2)ϕ(x + 2) dx, \\ N_{2} & = ∫_{0}^{∞} ϕ(x + 2)ϕ(x + 1) dx, \\ N_{3} & = ∫_{0}^{∞} ϕ(x + 1)ϕ(x + 1) dx. \end{aligned}

Numerical methods that compute N_{∗} in high precision are described in [11].

The mass matrix Y is similarly defined for {Y_{j}^{m}}, and the entries are

\begin{aligned}Y_{00} & = 3 − N_{1} − 2N_{2} − N_{3}, \\ Y_{01} & = y_{jr_{1}}^{1} + y_{jr_{2}}^{1}(1 − N_{1} − N_{2}) + y_{jr_{3}}^{1}(1 − N_{2} − N_{3}), \\ Y_{11} & = (y_{jr_{1}}^{1})^{2} + (y_{jr_{2}}^{1})^{2}(1 − N_{1}) + (y_{jr_{3}}^{1})^{2}(1 − N_{3}) − 2y_{jr_{2}}^{1}y_{jr_{3}}^{1}N_{2}. \end{aligned}

Here, for convenience, we define r_{0} = 2^{j} − 4, which is the last “interior” basis function at level j, and r_{i} = r_{0} + i, for i = 1, 2, 3. The value j will be clear in the context.

There exist S = X^{−1/2} and T = Y^{−1/2}. Let

( \begin{matrix}ϕ_{j,−1}(x) \\ ϕ_{j,0}(x) \end{matrix} ) = S ( \begin{matrix}X_{j}^{0}(x) \\ X_{j}^{1}(x) \end{matrix} ) , ( \begin{matrix}ϕ_{jr_{1}}(x) \\ ϕ_{jr_{2}}(x) \end{matrix} ) = T ( \begin{matrix}Y_{j}^{0}(x) \\ Y_{j}^{1}(x) \end{matrix} ) .

In the following we will denote the original compactly supported scaling functions by ϕ and the new scaling functions on [0, 1] by ϕ. One should understand that ϕ_{jk} = ϕ_{jk} for 1 ≤ k ≤ r_{0}. Then

{ϕ_{j,−1}, ϕ_{j0}, ϕ_{j1}, …, ϕ_{jr_{0}}, ϕ_{jr_{1}}, ϕ_{jr_{2}}}

is an orthonormal basis of V_{j}[0, 1]. See Figure 2.

Theoretically V_{j}[0, 1] preserve the nested relation [3]:

V_{j}[0, 1] ⊂ V_{j+1}[0, 1]

Figure 1. The nonzero part of X_{3}^{0}(x) (top left), X_{3}^{1}(x) (top right), Y_{3}^{0}(x) (bottom left) and Y_{3}^{1}(x).

Figure 2. The nonzero part of ϕ_{3,−1}(x) (top left), ϕ_{30}(x) (top right), ϕ_{3r_{1}}(x) (bottom left) and ϕ_{3r_{2}}(x).

However, there will be different sets of scaling equations for the “left,” “interior” and “right” basis functions. Those for the “interior” scaling functions are identically inherit from the original scaling equation (1.1). Those for the two edges will be numerically determined as following.

We now build the low pass filter, H. Let us illustrate the idea by the projection from V_{4}[0, 1] onto V_{3}[0, 1]. Let

Φ_{3} = ( ϕ_{3,−1} ; ϕ_{3,0} ; ϕ_{3,1} ; ⋮ ; ϕ_{3,4} ; ϕ_{3,5} ; ϕ_{3,6} ) , Φ_{4} = ( ϕ_{4,−1} ; ϕ_{4,0} ; ϕ_{4,1} ; ⋮ ; ϕ_{4,12} ; ϕ_{4,13} ; ϕ_{4,14} ) .

We will determine H such that

Φ_{3} = HΦ_{4}.

Now H = (H_{mn}) is an 8 × 16 matrix determined by

H_{mn} = (ϕ_{3m}, ϕ_{4n}).

We partition H into 9 parts according to the index m, n being in the “left,” “interior” or “right” portion of the basis functions.

H = ( \begin{matrix}H_{LL} & H_{LI} & H_{LR} \\ H_{IL} & H_{II} & H_{IR} \\ H_{RL} & H_{RI} & H_{RR} \end{matrix} )

We find H_{LL} is a 2 × 2 matrix with

\begin{aligned}H_{LL} & = ( \begin{matrix}(ϕ_{3,−1}, ϕ_{4,−1}) & (ϕ_{3,−1}, ϕ_{4,0}) \\ (ϕ_{3,0}, ϕ_{4,−1}) & (ϕ_{3,0}, ϕ_{4,0}) \end{matrix} ) \\ & = S ( \begin{matrix}(X_{3}^{0}, X_{4}^{0}) & (X_{3}^{0}, X_{4}^{1}) \\ (X_{3}^{1}, X_{4}^{0}) & (X_{3}^{1}, X_{4}^{1}) \end{matrix} ) S^{t}. \end{aligned}

H_{LI} is a 2 × 12 matrix which are zeros except for the first 3 columns:

\begin{matrix}( & \begin{matrix}(ϕ_{3,−1}, ϕ_{4,1}) & (ϕ_{3,−1}, ϕ_{4,2}) & (ϕ_{3,−1}, ϕ_{4,3}) \\ (ϕ_{3,0}, ϕ_{4,1}) & (ϕ_{3,0}, ϕ_{4,2}) & (ϕ_{3,0}, ϕ_{4,3}) \end{matrix} & ) \\ = & \frac{1}{\sqrt{2}} S ( \begin{matrix}c_{1} + c_{3} & c_{2} & c_{3} \\ x_{−1}^{1}c_{1} + x_{0}^{1}c_{3} & x_{0}^{1}c_{2} & x_{0}^{1}c_{3} \end{matrix} ) \end{matrix}

H_{II} is a familiar matrix of dimension 4 × 12.

\begin{aligned}H_{II} & = ( \begin{matrix}(ϕ_{3,1}, ϕ_{4,1}) & (ϕ_{3,1}, ϕ_{4,2}) & … & (ϕ_{3,1}, ϕ_{4,12}) \\ (ϕ_{3,2}, ϕ_{4,1}) & (ϕ_{3,2}, ϕ_{4,2}) & … & (ϕ_{3,2}, ϕ_{4,12}) \\ (ϕ_{3,3}, ϕ_{4,1}) & (ϕ_{3,3}, ϕ_{4,2}) & … & (ϕ_{3,3}, ϕ_{4,12}) \\ (ϕ_{3,4}, ϕ_{4,1}) & (ϕ_{3,4}, ϕ_{4,2}) & … & (ϕ_{3,4}, ϕ_{4,12}) \end{matrix} ) \\ & = \frac{1}{\sqrt{2}} ( \begin{matrix}0 & c_{0} & c_{1} & c_{2} & c_{3} & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & c_{0} & c_{1} & c_{2} & c_{3} & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & c_{0} & c_{1} & c_{2} & c_{3} & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & c_{0} & c_{1} & c_{2} & c_{3} & 0 \end{matrix} ) \end{aligned}

Similar to H_{LI}, H_{RI} are zeros except for the last three columns:

\begin{matrix}( & \begin{matrix}(ϕ_{3,5}, ϕ_{4,10}) & (ϕ_{3,5}, ϕ_{4,11}) & (ϕ_{3,5}, ϕ_{4,12}) \\ (ϕ_{3,6}, ϕ_{4,10}) & (ϕ_{3,6}, ϕ_{4,11}) & (ϕ_{3,6}, ϕ_{4,12}) \end{matrix} & ) \\ = \frac{1}{\sqrt{2}} T ( \begin{matrix}c_{0} & c_{1} & c_{2} + c_{0} \\ y_{5}^{1}c_{0} & y_{5}^{1}c_{1} & y_{5}^{1}c_{2} + y_{6}^{1}c_{0} \end{matrix} ) \end{matrix}

And H_{RR} is a 2 × 2 matrix with

\begin{aligned}H_{RR} & = ( \begin{matrix}(ϕ_{3,5}, ϕ_{4,13}) & (ϕ_{3,5}, ϕ_{4,14}) \\ (ϕ_{3,6}, ϕ_{4,13}) & (ϕ_{3,6}, ϕ_{4,14}) \end{matrix} ) \\ & = T ( \begin{matrix}(Y_{3}^{0}, Y_{4}^{0}) & (Y_{3}^{0}, Y_{4}^{1}) \\ (Y_{3}^{1}, Y_{4}^{0}) & (Y_{3}^{1}, Y_{4}^{1}) \end{matrix} ) T^{t} \end{aligned} .

By orthogonality, H_{LR}, H_{IL}, H_{IR}, H_{RL} are zero matrices.

In general, note that the dimension of H depends on the level index j; that is, 2^{j} × 2^{j+1}. However, the nonzero entries in the “boarder” parts of H do not depend on j, only the “interior” part H_{II} does. H_{LL} and H_{RR} are p × p matrices, and the nonzero parts of H_{LI} and H_{RI} are p×(2p−1). The “interior” part inherits the pyramid algorithm for the original compactly supported low pass filters [7].

Note that V_{j}[0, 1] ≠ V_{j}|_{[0,1]}. However,

V_{3}[0, 1] ⊂ V_{4}[0, 1] ⊂ ⋯

and ∪V_{j}[0, 1] is dense in L^{2}([0, 1]). The choice of δ_{L} = δ_{R} = 1 makes dim V_{j}[0, 1] = 2^{j}. Hence there will be an equal dimension for V_{j}[0, 1] and W_{j}[0, 1]. Functions in V_{j}[0, 1] are actually polynomials in the intervals [0, δ_{L}/2^{j}] and [1 − δ_{R}/2^{j}, 1], so they are smooth (Figure 2).

3. Wavelets on the interval. As in the framework of multiresolution approximation of L^{2}(ℝ), we now decompose

V_{j+1}[0, 1] = V_{j}[0, 1] ⊕ W_{j}[0, 1].

That is, define W_{j}[0, 1] to be the subspace of V_{j+1}[0, 1] which is orthogonal to V_{j}[0, 1]. Clearly we have

dim W_{j}[0, 1] = 2^{j}.

The wavelets on interval are the orthonormal basis functions for W_{j}[0, 1].

The original compactly supported wavelets satisfy

ψ_{j,k} = \frac{1}{\sqrt{2}} ∑_{l} g_{l−2k}ϕ_{j+1,l}.

Since g_{l−2k} ≠ 0 only for (−2 = 2 − 2p) ≤ l − 2k ≤ 1, thus

\{ \begin{matrix}1 ≤ l ≤ 2^{j+1} − 4 \\ −2 > −2k \\ 1 < 2^{j+1} − 4 − 2k \end{matrix}
⇒ 2 ≤ k ≤ (2^{j} − 3 = 2^{j} − p − 1)
⇒ ψ_{jk}(x) ∈ V_{j+1}[0, 1], for 2 ≤ k ≤ 2^{j} − 3.

But we would like to shift our mother wavelet such that suppψ(x) = [0, 2p − 1], thus the index would be shifted such that

ψ_{jk}(x) ∈ V_{j+1}[0, 1], for 1 ≤ k ≤ 2^{j} − 4 = r_{0}.

That is, we already have 2^{j} − 4 orthonormal basis on W_{j}[0, 1], which are the original compactly supported wavelets. We need 4 = 2p more basis functions for W_{j}[0, 1]; they are selected such that there are 2 = p basis functions near each boundary point.

From above, we know the orthonormal bases for V_{j}[0, 1] and V_{j+1}[0, 1], and we have 2^{j} −4 basis functions on W_{j}[0, 1]. Now we want the basis functions of V_{j+1}[0, 1] which do not belong to the span of V_{j}[0, 1] and {ψ_{j,k}}_{k=1}^{r_{0}}. Then we find their components which are orthogonal to V_{j}[0, 1]. The result will be functions in W_{j}[0, 1].

We can show that

ϕ(x)_{j+1,−1} ∉ span{V_{j}[0, 1] ∪ {ψ_{jk}}_{k=1}^{r_{0}}},
ϕ(x)_{j+1,0} ∉ span{V_{j}[0, 1] ∪ {ψ_{jk}}_{k=1}^{r_{0}} ∪ ϕ(x)_{j+1,−1}}.

Hence, on the “left,” we choose the orthogonal complements of ϕ_{j+1,−1} and ϕ_{j+1,0}; similarly, on the “right,” we choose ϕ_{j+1,r_{1}} and ϕ_{j+1,r_{2}}. Define the complement functions:

ϕ_{j0}^{L} = ϕ_{j+1,−1} − P_{j}ϕ_{j+1,−1}
ϕ_{j1}^{L} = ϕ_{j+1,0} − P_{j}ϕ_{j+1,0}
ϕ_{j0}^{R} = ϕ_{j+1,r_{1}} − P_{j}ϕ_{j+1,r_{1}}
ϕ_{j1}^{R} = ϕ_{j+1,r_{2}} − P_{j}ϕ_{j+1,r_{2}}

Here P_{j} is the L^{2}-orthogonal projection from L^{2}([0, 1]) onto V_{j}[0, 1]. Then ϕ_{j∗}^{∗} ∈ W_{j}[0, 1] and it can be shown that they are linearly independent. They are already orthogonal with the “interior” wavelets {ψ_{jk}}_{k=1}^{r_{0}}. Now we orthonormalize ϕ_{j∗}^{∗} to get a set of orthonormal basis for W_{j}[0, 1].

As for the orthonormalization for scaling functions on the “edge,” we need the mass matrix of ϕ_{j∗}^{L}, say X^{⊥}, and that of ϕ_{j∗}^{R}, say Y^{⊥}. It turns out that we do not need much more calculations,

X^{⊥} = I − H_{LL}^{t}H_{LL}, Y^{⊥} = I − H_{RR}^{t}H_{RR}.

Here I is the identity matrix of compatible dimension, eg, 2 × 2 for this case. Again, X^{⊥} and Y^{⊥} are symmetric positive definite and they are independent of j. Now define

( \begin{matrix}ψ_{j,−1}(x) \\ ψ_{j,0}(x) \end{matrix} ) = (X^{⊥})^{− \frac{1}{2}} ( \begin{matrix}ϕ_{j0}^{L}(x) \\ ϕ_{j1}^{L}(x) \end{matrix} ) , ( \begin{matrix}ψ_{j,r_{1}}(x) \\ ψ_{j,r_{2}}(x) \end{matrix} ) = (Y^{⊥})^{− \frac{1}{2}} ( \begin{matrix}ϕ_{j0}^{R}(x) \\ ϕ_{j1}^{R}(x) \end{matrix} ) .

In the following we will denote the original compactly supported wavelets by ψ and the new wavelets on [0, 1] by ψ. One should understand that ψ_{jk} = ψ_{jk} for 1 ≤ k ≤ r_{0}. Then

{ψ_{j,−1}, ψ_{j0}, ψ_{j1}, …, ψ_{jr_{0}}, ψ_{jr_{1}}, ψ_{jr_{2}}}

forms an orthonormal basis of W_{j}[0, 1]. See Figure 3.

For the relation between V_{j+1}[0, 1] and W_{j}[0, 1], it amounts to construct the high pass filter, G. Again, we illustrate the idea by the projection from V_{4}[0, 1] onto W_{3}[0, 1]. Let

Ψ_{3} = ( ψ_{3,−1} ; ψ_{3,0} ; ψ_{3,1} ; ⋮ ; ψ_{3,4} ; ψ_{3,5} ; ψ_{3,6} ) , Φ_{4} = ( ϕ_{4,−1} ; ϕ_{4,0} ; ϕ_{4,1} ; ⋮ ; ϕ_{4,12} ; ϕ_{4,13} ; ϕ_{4,14} )

We will determine G such that

Ψ_{3} = GΦ_{4}.

Now G is a (8 = 2^{j}) × (16 = 2^{j+1}) matrix. We partition G in the same pattern as for H,

G = ( \begin{matrix}G_{LL} & G_{LI} & G_{LR} \\ G_{IL} & G_{II} & G_{IR} \\ G_{RL} & G_{RI} & G_{RR} \end{matrix} ) .

After some manipulations, we can write G in terms of H:

G_{LL} = (I − H_{LL}^{t}H_{LL})^{\frac{1}{2}}
G_{LI} = −(I − H_{LL}^{t}H_{LL})(H_{LL}^{t}H_{LI})
G_{RI} = −(I − H_{RR}^{t}H_{RR})(H_{RR}^{t}H_{RI})
G_{RR} = (I − H_{RR}^{t}H_{RR})^{\frac{1}{2}}

The remarks for H also apply to G. Hence all but a fixed small number (2p) of filter coefficients are the same as the original compactly supported wavelets. From the construction, we have

H^{t}H + G^{t}G = I.

Therefore, the fast pyramid decomposition and reconstruction algorithm of Mallat [7] is almost totally preserved.

4. Pre-filters. In the application of discrete wavelet transforms, the signals (or images) are in discrete values; namely in form of vectors or matrices. In this case, as pointed out in [4], the discrete wavelets transform on the interval has a shortcoming. The decomposed data may have larger high frequency components than the periodized discrete wavelet transform, as illustrated in Figures 4.

Let us explain the situation for the case of p = 2. For wavelets of this approximation order, polynomials of degree 0 and 1 are in V_{j}, in the following rough sense:

1 = ∑_{k}(1, ϕ_{jk})ϕ_{jk} x = ∑_{k}(x, ϕ_{jk})ϕ_{jk}

Figure 3. The nonzero part of ψ_{3,−1}(x) (top left), ψ_{30}(x) (top right), ψ_{3r_{1}}(x) (bottom left) and ψ_{3r_{2}}(x).

Figure 4. If we take 1024 samples from f(x) = x + cos x in [0, 1] and perform the periodized discrete wavelet transform (DWT), that is, extend f(x) periodically on ℝ and perform the DWT. There are artificial high amplitudes in each W_{j} due to the discontinuity introduced by the periodic extension (top). When the interval DWT is applied without the filtering process, the amplitudes are even larger (middle). However, the interval DWT with filtering process results flat high frequency signals. It is because cos x is smooth and the (processed) high pass filter annihilates x.

For wavelets on the interval, it is true that we still have:

1 = ∑_{k=−1}^{2^{j}−2}(1, ϕ_{jk})ϕ_{jk} x = ∑_{k=−1}^{2^{j}−2}(x, ϕ_{jk})ϕ_{jk}

However, (1, ϕ_{jk}) = 1 and (x, ϕ_{jk}) = k + M_{1}. Hence the vectors corresponding to the projection coefficients of 1 and x in V_{0} are essentially

(4.1)  x_{0} = (…, 1, 1, 1, …)^{t} and x_{1} = (…, 1, 2, 3, …)^{t}.

While these two vectors can be annihilated by the high pass filter of the discrete wavelet transform (DWT), it is no longer true for our new DWT based on wavelets on the interval. In practice, signals are pointwise values of a function whose most significant parts are in the polynomials of low degree; that is, in the form of

α_{0}x_{0} + α_{1}x_{1} + ⋯,

where α_{0} and α_{1} are relatively large. The periodized high pass filter annihilates entire x_{0} and most part of x_{1}, while the interval high pass filter fails to do that. Sometimes, as we see in Figure 4, it may introduce very large high frequency components (very large coefficients in W_{j}[0, 1] for high j).

However, the projection coefficients of 1 and x onto V_{j}[0, 1] are still in the form of (4.1), except for the first and last two (p) components. The idea of the prefilter in [4] is to transform the space span{x_{0}, x_{1}} into which can be annihilated by the high pass filter G. That is, essentially

( 1 1 1 ⋮ 1 1 1 ) + ( 1 ; 2 ; 3 ; ⋮ ; 2^{j} − 2 ; 2^{j} − 1 ; 2^{j} ) ↦ ( ∫ ϕ_{j,−1} ; ∫ ϕ_{j0} ; ∫ ϕ_{j1} ; ⋮ ; ∫ ϕ_{j,r_{0}} ; ∫ ϕ_{j,r_{1}} ; ∫ ϕ_{j,r_{2}} ) + ( ∫ xϕ_{j,−1} ; ∫ xϕ_{j0} ; ∫ xϕ_{j1} ; ⋮ ; ∫ xϕ_{j,r_{0}} ; ∫ xϕ_{j,r_{1}} ; ∫ xϕ_{j,r_{2}} )

Since the spaces are different only at the first and last two entries, this amounts to find linear operators F_{L} and F_{R} such that

F_{L} : span\{( \begin{matrix}∫ ϕ_{j,−1} \\ ∫ ϕ_{j,0} \end{matrix} ) , ( \begin{matrix}∫ xϕ_{j,−1} \\ ∫ xϕ_{j,0} \end{matrix} )\} ↦ span\{( \begin{matrix}∫ ϕ_{j,−1} \\ ∫ ϕ_{j,0} \end{matrix} ) , ( \begin{matrix}∫ xϕ_{j,−1} \\ ∫ xϕ_{j,0} \end{matrix} )\}

F_{R} : span\{( \begin{matrix}∫ ϕ_{j,2^{j}−3} \\ ∫ ϕ_{j,2^{j}−2} \end{matrix} ) , ( \begin{matrix}∫ xϕ_{j,2^{j}−3} \\ ∫ xϕ_{j,2^{j}−2} \end{matrix} )\} ↦ span\{( \begin{matrix}∫ ϕ_{j,2^{j}−3} \\ ∫ ϕ_{j,2^{j}−2} \end{matrix} ) , ( \begin{matrix}∫ xϕ_{j,2^{j}−3} \\ ∫ xϕ_{j,2^{j}−2} \end{matrix} )\}

Since the vectors listed above are linearly independent, so F_{L} and F_{R} will be non-singular matrices. After some linear combinations of the basis functions in these subspaces, we observe that F_{L} and F_{R} are actually independent of j:

\begin{aligned}F_{L} & = SX̃ ( \begin{matrix}K_{2}^{0} & −M_{1}K_{2}^{0} + K_{2}^{1} \\ K_{1}^{0} & −M_{1}K_{1}^{0} + K_{1}^{1} \\ 1 & 0 \end{matrix} ) ( \begin{matrix}1 & −1 \\ 1 & 0 \end{matrix} )^{−1} \\ F_{R} & = TỸ ( \begin{matrix}1 & −1 \\ 1 − K_{2}^{0} & K_{2}^{0}(M_{1} − 2) − K_{2}^{1} \\ 1 − K_{1}^{0} & K_{1}^{0}(M_{1} − 2) − K_{1}^{1} + 1 \end{matrix} ) ( \begin{matrix}1 & −1 \\ 1 & 0 \end{matrix} )^{−1} \end{aligned}

where

X̃ = ( \begin{matrix}x_{−2}^{0} & x_{−1}^{0} & x_{0}^{0} \\ x_{−2}^{1} & x_{−1}^{1} & x_{0}^{1} \end{matrix} ) , Ỹ = ( \begin{matrix}y_{2^{j}−2}^{0} & y_{2^{j}−2}^{0} & y_{2^{j}−2}^{0} \\ y_{2^{j}−2}^{1} & y_{2^{j}−1}^{1} & y_{2^{j}−1}^{1} \end{matrix} ) .

Here K_{n}^{m} = ∫_{0}^{∞} x^{m}ϕ(x + n) dx. Note that m ≥ 0 and 1 ≤ n ≤ (2 = 2p − 2) are the meaningful indices. We also follow [11] for high accurate values of K_{n}^{m}.

Now we have the 2×2 (or p×p) filter matrices F_{L} and F_{R}. For one-dimensional data (signals), we preprocess the input v = (v_{0}, v_{1}, …, v_{n})^{t} by converting it to ṽ = (ṽ_{0}, ṽ_{1}, v_{2}, …, v_{n−2}, ṽ_{n−1}, ṽ_{n})^{t}, where

( ṽ_{0} ; ṽ_{1} ) = F_{L} ( \begin{matrix}v_{0} \\ v_{1} \end{matrix} ) , ( ṽ_{n−1} ; ṽ_{n} ) = F_{R} ( \begin{matrix}v_{n−1} \\ v_{n} \end{matrix} ) .

Then we perform the interval DWT on ṽ. Reversely, we post-process ṽ by F_{L}^{−1} and F_{R}^{−1} after the inverse transform. See also Figure 4 for a demonstration. There is almost no operation count added by the pre- and post-process.

For two-dimensional data (images), however, we need the preprocess on the two “boarder” rows and columns. Assume that one level of the two-dimensional wavelet transform is implemented by applying one level of the one-dimensional wavelet transform on each row of data, followed by one level of the one-dimensional wavelet transform on each column of data. See [7]. In this case, multiplying F_{L} (F_{R}) to the first (last) two rows and columns will not work for our purpose.

The key point of the filtering process is to keep the polynomials 1, x, y and xy invariant under the low pass filter, while orthogonal to the high pass filter. Keep this in mind, we find a successful filtering process. That is, we apply the prefilters to column vectors, perform interval DWT on each row, undo the prefilters (switch it back to the original space), then apply the prefilters to row vectors and perform the interval DWT on each column, then undo the prefilters. The inverse transform needs reversed processes. For coding details, we refer to the master dissertation of the first author, Miss Cheng [2]. Figure 5 shows the effect of this filtering process.

Our filtering process introduces O(pN) of operation count to the computational complexity, where N is the dimension of input matrix. Whether it worths the efforts is yet to be seen by more experiments on actual image samples. More can be found in [2].

REFERENCES

[1] L. Andersson, N. Hall, B. Jawerth and G. Peters, Wavelets on closed subsets of the real line, in “Recent advances in wavelet analysis,” Academic Press, Boston, 1994, pp. 1–61.

[2] J. Cheng, The construction of wavelets on the interval and a prefilter process for wavelet transforms, Master Dissertation, Department of Mathematics, National Central University, 1995.

[3] A. Cohen, I. Daubechies, B. Jawerth and P. Vial, Multiresolution analysis, wavelets and fast algorithms on the interval, C. R. Acad. Sci. Paris, Série I, 316 (1992), pp. 417–421.

[4] A. Cohen, I. Daubechies and P. Vial, Wavelets on the interval and fast wavelet transforms, Appl. Comput. Harmonic Anal., 1 (1993), pp. 54–81.

[5] I. Daubechies, Orthonormal bases of compactly supported wavelets, Comm. Pure Appl. Math., 41 (1988), pp. 909–996.

[6] I. Daubechies, “Ten lectures on wavelets,” SIAM, Philadelphia, 1992.

[7] S. Mallat, Multiresolution approximations and wavelet orthonormal bases of L^{2}(ℝ), Trans. Amer. Math. Soc., 315 (1989), pp. 69–87.

[8] S. Mallat, A theory for multiresolution signal decomposition: The wavelet representation, IEEE Trans. on Patt. Anal. Mach. Intell., 11 (1989), pp. 674–693.

[9] Y. Meyer, “Ondelletes et fonctions splines,” Seminaire Equations aux Derivees Partielles, Ecole Polytechnique, Paris, 1986.

[10] Y. Meyer, Ondelettes sur l’intervalle, Rev. Math. Iberoamer., 7 (1992), pp. 115–133.

[11] W.-C. Shann and C.-C. Yen, Quadratures involving polynomials and Daubechies’ wavelets, Technical Report 9301, Department of Mathematics, National Central University, 1993.

[12] G. Strang, Wavelets and dilation equations: A brief introduction, SIAM Review, 31 (1989), pp. 614–627.

Figure 5. We choose the test two-dimensional data to be the 256 × 256 matrix A = (a_{ij}) where a_{ij} = −1 if i = j, and a_{ij} = \frac{1}{\sqrt{|i−j|}} if i ≠ j. Take periodized 2-d DWT of A, we show only

the entries whose magnitude is larger than 0.1% of the maximum (top). Observe the black edges caused by the discontinuity introduced by the periodic extension. Similar to the one-dimensional case, the interval 2-d DWT has worse result (bottom left). The interval 2-d DWT with the filtering process completely eliminates the artificial large entries (bottom right). If we set the entries not shown on the graph by 0, then the compression rate has raised from 10.7 (top) to 14.3 (bottom right).