# Introduction to wavelets in engineering

John Williams, Kevin Amaratunga

► To cite this version:

John Williams, Kevin Amaratunga. Introduction to wavelets in engineering. International Journal for Numerical Methods in Engineering, 1994, 37, pp.2365-2388. 10.1002/nme.1620371403. hal-01311772

HAL Id: hal-01311772

https://hal.science/hal-01311772

Submitted on 4 May 2016

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>
# INTRODUCTION TO WAVELETS IN ENGINEERING

JOHN R. WILLIAMS* AND KEVIN AMARATUNGA†

*Intelligent Engineering Systems Laboratory, Massachusetts Institute of Technology, Cambridge, MA 02139, U.S.A.*

## SUMMARY

The aim of this paper is to provide an introduction to the subject of wavelet analysis for engineering applications. The paper selects from the recent mathematical literature on wavelets the results necessary to develop wavelet-based numerical algorithms. In particular, we provide extensive details of the derivation of Mallat's transform and Daubechies' wavelet coefficients, since these are fundamental to gaining an insight into the properties of wavelets. The potential benefits of using wavelets are highlighted by presenting results of our research in one- and two-dimensional data analysis and in wavelet solutions of partial differential equations.

## 1. INTRODUCTION

### 1.1. *History of wavelets*

The appearance of wavelets is a relatively recent development in mathematics. The name wavelet or *ondelette* was coined some ten years ago by French researchers, including Morlet *et al.*,<fnref n="1" /> Morlet² and Grossmann and Morlet.³ The existence of wavelet-like functions has been known since the early part of the century (notable examples being what are now known as the Haar wavelet and the Littlewood-Paley wavelet) and many of the ideas now embodied in wavelets originated from work in subband coding in engineering, coherent states and renormalization group theory in physics and the study of Calderon-Zygmund operators in mathematics. However, it was not until recently that the unifying concepts necessary for a general understanding of wavelets were provided by researchers such as Stromberg,⁴ Meyer,⁵ Mallat,⁶ Daubechies,⁷ Battle⁸ and Lemarié.⁹ Since then, the growth of wavelet research in mathematics has been explosive with numerous authors contributing significantly.

In this paper we attempt to provide an introduction to wavelets which is attuned to the needs of the engineering researcher who may wish to apply wavelet techniques to his or her domain. We have found that the biggest *barrier to entry* is the derivation of the equations needed to program the algorithms. We have therefore provided in great detail the derivations of the matrices used in Mallat's wavelet transform, in the development of the wavelet coefficients and in the solution of partial differential equations.

### 1.2. *Some properties of wavelets*

The wavelet transform is a tool that cuts up data, functions or operators into different frequency components, and then studies each component with a resolution matched to its scale.*

<footnote>* Associate Professor</footnote>

<footnote>† Graduate Student</footnote>

<footnote n="1">¹ Daubechies¹⁰</footnote>
For example, in signal analysis, the wavelet transform allows us to view a time history in terms of its frequency components. In this respect, the wavelet transform performs in much the same way as the Fourier transform. The Fourier transform extracts from the signal details of the frequency content but loses all information on the location of a particular frequency within the signal. Time localization must then be achieved by first windowing the signal, and then by taking its Fourier transform. The problem with windowing is that the slice of the signal that is extracted is always the same length. Thus, the time slice (number of data points) used to resolve a high-frequency component is the same as the number used to resolve a low-frequency component.

In contrast to windowed Fourier transforms, the wavelet adapts the width of its time slice according to the frequency components being extracted. For example, if we express a windowed Fourier transform as

$$ (T^{\text{win}}f)(\omega, t) = \int f(s)g(s-t)e^{-i\omega t} ds \quad (1) $$

then the corresponding wavelet transform is

$$ (T^{\text{wav}}f)(a, b) = \int f(t)\psi\left(\frac{t-b}{a}\right)dt \quad (2) $$

where the position of the slice of the signal that the wavelet samples in time is controlled by *b* and the extent of the slice by the scaling parameter *a*.

Wavelets are being used in signal processing to develop *best basis* methods for signal compression and for signal characterization.

Wavelet bases can also be used to precondition stiffness matrices and to develop fast hierarchical solution procedures for partial differential equations.

In Section 2 we describe the general properties that we require of scaling functions and wavelets, using the Haar system as a specific example.

In Section 3 we show how to design and construct Daubechies' orthogonal compactly supported wavelet functions. The properties of other wavelet families are outlined in Section 4.

Section 5 outlines the theory of wavelet multiresolution analysis developed by Mallat.⁶,¹¹

In Sections 6 and 7 we show how wavelets may be used to solve engineering problems including: (1) analysis of image data and (2) solution of partial differential equations.

## 2. MULTILEVEL REPRESENTATION OF A FUNCTION

The goal of multiresolution analysis is to develop representations of a function $f(x)$ at various levels of resolution. To achieve this we seek to expand the given function in terms of basis functions $\phi(x)$ which can be scaled to give multiple resolutions of the original function.

In order to develop a multilevel representation of a function in $L^2(\mathbf{R})$ we seek a sequence of embedded subspaces $V_i$ such that

$$ \{0\} \cdots \subset V_{-1} \subset V_0 \subset V_1 \subset V_2 \cdots \subset L^2(\mathbf{R}) \quad (3) $$

with the following properties:

(a) $\bigcup_{j \in \mathbb{Z}} V_j$ is dense in $L^2(\mathbf{R})$. (We use $\mathbb{Z}$ to denote the set of integers.)

(b) $\bigcap_{j \in \mathbb{Z}} V_j = \{0\}$.

(c) The embedded subspaces are related by a scaling law

$$ g(x) \in V_j \Leftrightarrow g(2x) \in V_{j+1} $$
(d) Each subspace is spanned by integer translates of a single function $g(x)$ such that

$$g(x) \in V_0 \Leftrightarrow g(x+1) \in V_0$$

We then seek a scaling function $\phi(x) \in V_0$ such that its integer translates $\{\phi(x-k), k \in \mathbb{Z}\}$ from a Riesz basis for the space $V_0$. Then $\phi(2x-k)$ form a basis for the space $V_1$. (We note that scaling by powers other than two is possible but non-powers of 2 will not be discussed here.) Thus,

$$\begin{aligned} V_0 &= \operatorname{span} \frac{\{\phi(x-k), k \in \mathbb{Z}\}}{} \\ V_1 &= \operatorname{span} \frac{\{\phi(2x-k), k \in \mathbb{Z}\}}{} \end{aligned}$$

Since the space $V_0$ lies within the space $V_1$, we can express any function in $V_0$ in terms of the basis functions of $V_1$. In particular,

$$\phi(x) = \sum_{k=-\infty}^{\infty} a_k \phi(2x-k) \quad (4)$$

in which $a_k, k \in \mathbb{Z}$ is a square summable sequence. We will refer to equation (4) as the *dilation equation*. The dilation equation is sometimes referred to as the *refinement equation* or the *scaling relation*.

If we now define

$$\phi_{m,k}(x) = 2^{m/2} \phi(2^m x - k) \quad (5)$$

then $\phi_{m,k}(x)$, $k \in \mathbb{Z}$ forms a Riesz basis for the space $V_m$. We shall refer to the dilation parameter $m$ as the *scale*.

## 2.1. Example using the box function as a basis

Let us demonstrate the validity of our assumptions so far by using the box function (Figure 1(a)) as our basis function $\phi(x)$,

$$\phi(x) = \begin{cases} 1, & 0 \le x < 1 \\ 0, & \text{otherwise} \end{cases} \quad (6)$$

From Figure 1(b), it is evident that $\{\phi(x-k), k \in \mathbb{Z}\}$ is an orthonormal set. We may now (in a crude sense) approximate a function $f \in L^2(\mathbb{R})$ by its projection $P_0f$ onto the space $V_0$:

$$P_0f = \sum_{k=-\infty}^{\infty} c_{0,k} \phi(x-k) \quad (7)$$

In general, a function may be approximated by its projection $P_mf$ onto the space $V_m$:

$$P_m f = \sum_{k=-\infty}^{\infty} c_{m,k} \phi_{m,k}(x) \quad (8)$$

and in fact $P_m f$ approaches $f$ as $m \to \infty$. Figure 2(b) and 2(d) shows two such approximations to a function at consecutive scales.

We note that the box function satisfies the dilation equation with coefficients $a_0 = a_1 = 1$,

$$\phi(x) = \phi(2x) + \phi(2x-1) \quad (9)$$

from which we see that

$$V_0 \subset V_1 \subset V_2 \subset \dots$$
<figure><img src="image_2.png" /><figcaption>Figure 1. (a) The Haar scaling function $\phi(x) = \phi_{0,0}$. (b) Translations of the Haar scaling function $\phi_{0,0}$, $\phi_{0,1}$ and $\phi_{0,2}$. (c) Dilations of the Haar scaling function $\phi_{-1,0}$ and $\phi_{-1,1}$. (d) The Haar wavelet $\psi(x) = \psi_{0,0}$.</figcaption></figure>

<figure><img src="image_3.png" /><figcaption>Figure 2. Successive approximations to a function. (a) $f = \sin(2\pi x) \sin(4\pi x)$. (b) Approximation at scale 4, $P_4 f$. (c) Detail at scale 4, $Q_4 f$. (d) Approximation at scale 5, $P_5 f$.</figcaption></figure>
In this example we have shown that the box function can be used to develop a sequence of embedded spaces with each space spanned by the translates of the box function at that scale.

## 2.2. Wavelets

Now let us investigate the difference between subspaces $V_{m-1}$ and $V_m$. Let us define a new subspace $W_{m-1}$ such that it is the orthogonal complement of $V_{m-1}$ in $V_m$,

$$V_m = V_{m-1} \oplus W_{m-1}, \quad V_{m-1} \perp W_{m-1} \qquad (10)$$

where $\oplus$ represents a direct sum. It follows then that the spaces $W_j$ are orthogonal and that

$$\bigoplus_{j \in \mathbb{Z}} W_j = \mathbf{L}^2(\mathbf{R}) \qquad (11)$$

Now let us introduce a wavelet function $\psi(x)$ such that $\psi(x-k)$ form a Riesz basis for the subspace $W_0$. Then

$$\psi_{m,k} = 2^{m/2} \psi(2^m x - k)$$

is a Riesz basis for $W_m$. If, in addition, the $\{\psi(x-k), k \in \mathbb{Z}\}$ form an orthonormal set, then it follows that $\{\psi_{m,k}, m, k \in \mathbb{Z}\}$ form an orthonormal basis for $\mathbf{L}^2(\mathbf{R})$.

Let us denote the projection of $f$ on $W_m$ as $Q_m f$. Then, from equation (10), we have

$$P_m f = P_{m-1} f + Q_{m-1} f \qquad (12)$$

This means that $Q_m f$ represents the detail that needs to be added to get from one level of approximation to the next finer level of approximation.

Furthermore, since the space $W_0$ is contained in the space $V_1$, we can express the wavelet function in terms of the scaling function at the next higher scale,

$$\psi(x) = \sum_{k=-\infty}^{\infty} b_k \phi(2^k x - k) \qquad (13)$$

The wavelet which corresponds to the box function is called the Haar wavelet and is shown in Figure 1(d). This wavelet satisfies equation (13) with coefficients $b_0 = 1$ and $b_1 = -1$, so that if $\phi(x)$ is the box function,

$$\psi(x) = \phi(2x) - \phi(2x-1). \qquad (14)$$

## 2.3. Multiresolution analysis

Multiresolution decomposition takes the expansion coefficients $c_{m,k}$ of an approximation, $P_m f$, to a function at scale $m$ and decomposes them into

(1) the expansion coefficients, $c_{m-1,k}$, of the approximation $P_{m-1} f$ at the next coarser scale $m-1$,

(2) the expansion coefficients, $d_{m-1,k}$, of the detail component $Q_{m-1} f = P_m f - P_{m-1} f$.

Section 5 details the algorithm for this decomposition. We can then repeat this process on the expansion coefficients $c_{m-1,k}$ of the approximation $P_{m-1} f$ to find the coefficients $c_{m-2,k}$ and $d_{m-2,k}$, and so on. Multiresolution decomposition therefore breaks down the original $\mathbf{L}^2$ space into a series of orthogonal subspaces at different resolutions.
### 3. CONSTRUCTION OF WAVELET SYSTEMS

Here we consider the problem of generating other wavelet functions with the special properties of orthogonality that we have seen in the case of the Haar function. Historically, special cases of wavelet functions, and functions with some of the properties of wavelets, were known long before Daubechies, Mallat and Meyer provided a general theory of wavelet functions.

#### 3.1. Generalized scaling functions and wavelets

The Haar scaling function and the Haar wavelet described in Section 2 are special cases of a more general class of functions. In general, a scaling function, $\phi(x)$, is the solution to a dilation equation of the form

$$ \phi(x) = \sum_{k=-\infty}^{\infty} a_k \phi(Sx - k) \quad (15) $$

A convenient choice of the dilation factor is $S = 2$, in which case the equation becomes

$$ \phi(x) = \sum_{k=-\infty}^{\infty} a_k \phi(2x - k) \quad (16) $$

The constant coefficients $a_k$ are called filter coefficients and it is often the case that only a finite number of these are non-zero. The filter coefficients are derived by imposing certain conditions on the scaling function. One of these conditions is that the scaling function and its translates should form an orthonormal set, i.e.

$$ \int_{-\infty}^{\infty} \phi(x) \phi(x+l) \, dx = \delta_{0,l}, \quad l \in \mathbb{Z} $$

where

$$ \delta_{0,l} = \begin{cases} 1, & l=0 \\ 0, & \text{otherwise} \end{cases} $$

A wavelet, $\psi(x)$, is orthogonal to the scaling function and is defined by

$$ \psi(x) = \sum_{k=-\infty}^{\infty} (-1)^k a_{N-1-k} \phi(2x-k) \quad (17) $$

where $N$ is an even integer.† This definition satisfies orthogonality, since

$$ \begin{aligned} \langle \phi(x), \psi(x) \rangle &= \int_{-\infty}^{\infty} \sum_{k=-\infty}^{\infty} a_k \phi(2x-k) \sum_{l=-\infty}^{\infty} (-1)^l a_{N-1-l} \phi(2x-l) \, dx \\ &= \frac{1}{2} \sum_{k=-\infty}^{\infty} (-1)^k a_k a_{N-1-k} \\ &= 0. \end{aligned} $$

The sets of coefficients $\{a_k\}$ and $\{(-1)^k a_{N-1-k}\}$ are said to form a pair of *quadrature mirror filters*.

<footnote>† Some texts define the filter coefficients of the wavelet as $(-1)^k a_{1-k}$. Equation (17), however, is more convenient to use when there are only a finite number of filter coefficients $a_0, \dots, a_{N-1}$, since it leads to a wavelet that has support over the same interval, $[0, N-1]$, as the scaling function.</footnote>
Figure 3 illustrates four scaling functions and their corresponding wavelets. These have been obtained by using different sets of filter coefficients to solve equations (16) and (17).

## 3.2. Derivation of filter coefficients

If the scaling function is to form a useful basis for functional analysis, it should display certain properties which lead to corresponding conditions on the filter coefficients. Consider the following properties:

(i) In order to uniquely define all scaling functions of a given shape, the area under the scaling function is normalized to unity, i.e.

$$ \int_{-\infty}^{\infty} \phi(x) dx = 1 \qquad (18) $$

The fact that the scaling function has a non-vanishing integral leads to the following condition on the filter coefficients:

$$ \sum_{k=-\infty}^{\infty} a_k = 2 \qquad (19) $$

(ii) For the scaling function to be orthogonal to its integer translates, the filter coefficients must satisfy the additional requirement that

$$ \int_{-\infty}^{\infty} \phi(x) \phi(x+l) dx = \delta_{0,l}, \quad l \in \mathbb{Z} \qquad (20) $$

This yields the condition

$$ \sum_{k=-\infty}^{\infty} a_k a_{k+2l} = 2\delta_{0,l}, \quad l \in \mathbb{Z} \qquad (21) $$

Once again, equation (21) is only a necessary condition for equation (20) to be satisfied. However, the condition is almost sufficient to ensure orthogonality and in practice it may be regarded as being sufficient.

(iii) Equations (19) and (21) are insufficient to determine a unique set of filter coefficients. In an N coefficient system, they yield a total of N/2 + 1 equations. Another N/2 - 1 equations are therefore required for a unique solution. One way to obtain the remaining equations, which will be considered here, is to require the scaling function to be able to exactly represent polynomials of order up to, but not greater than, p. Enforcing this requirement leads to the compactly supported wavelets developed by Daubechies.⁷ It will be seen shortly that p = N/2. The requirement for approximation of order p is that any function of the form

$$ f(x) = \alpha_0 + \alpha_1 x + \alpha_2 x^2 + \cdots + \alpha_{p-1} x^{p-1} \qquad (22) $$

can be exactly represented by an expansion of the form

$$ f(x) = \sum_{k=-\infty}^{\infty} c_k \phi(x-k) \qquad (23) $$

This may be translated into a condition on the wavelet. Taking the inner product of equation (23) with $\psi(x)$ gives

$$ \langle f(x), \psi(x) \rangle = \sum_{k=-\infty}^{\infty} c_k \langle \phi(x-k), \psi(x) \rangle = 0 $$
<figure><img src="image_4.png" /><figcaption>Figure 3. (a) Haar scaling function, (b) Haar wavelet, (c) Daubechies D4 scaling function, (d) Daubechies D4 wavelet, (e) Daubechies D6 scaling function, (f) Daubechies D6 wavelet, (g) Battle-Lemarié scaling function, (h) Battle-Lemarié wavelet</figcaption></figure>
Thus, from equation (22),

$$ \alpha_0 \int_{-\infty}^{\infty} \psi(x) dx + \alpha_1 \int_{-\infty}^{\infty} \psi(x) x dx + \cdots + \alpha_{p-1} \int_{-\infty}^{\infty} \psi(x) x^{p-1} dx \equiv 0 $$

This identity is valid for all $\alpha_j (j = 0, 1, 2, \dots, p-1)$. Choosing $\alpha_l = 1$ and all other $\alpha_j = 0$ gives

$$ \int_{-\infty}^{\infty} \psi(x) x^l dx = 0, \quad l = 0, 1, 2, \dots, p-1 \qquad (24) $$

Thus the first $p$ moments of the wavelet must be zero. To determine what constraints this places on the filter coefficients, substitute equation (17) in equation (24). Some calculation shows that

$$ \sum_{k=-\infty}^{\infty} (-1)^k a_k k^l = 0, \quad l = 0, 1, 2, \dots, p-1 \qquad (25) $$

It has already been shown that equation (25) must yield $N/2 - 1$ new equations. The equation obtained by putting $l = 0$ is redundant since it may be obtained from

$$ \sqrt{2[\text{equation (21)}]_{l=0} + 4[\text{equation (21)}]_{l>0} - [\text{equation (19)}]^2} $$

Thus, equation (25) yields $p-1$ new equations which means that $p = N/2$. Then the final condition on the filter coefficients becomes

$$ \sum_{k=-\infty}^{\infty} (-1)^k a_k k^l = 0, \quad l = 0, 1, 2, \dots, N/2 - 1 \qquad (26) $$

The filter coefficients $\{a_k, k = 1, \dots, N-1\}$ for an $N$ coefficient system are uniquely defined by equations (19), (21) and (26). Table I lists the coefficients for the Daubechies D4, D6 and D8 wavelet systems.

The problem which now remains is to construct the actual scaling function $\phi(x)$ from its filter coefficients.

### 3.3. Construction of scaling functions

In general, scaling functions do not have a closed-form solution. Instead, they have to be generated recursively from the dilation equation (equation (16)). The dilation equation may be

<table><caption>Table I. Daubechies filter coefficients for $N = 4, 6$ and $8$</caption><thead><tr><th>k</th><th>D4</th><th>D6</th><th>D8</th></tr></thead><tbody><tr><td>0</td><td>0·68301270189222</td><td>0·47046720778405</td><td>0·32580342805130</td></tr><tr><td>1</td><td>1·18301270189222</td><td>1·14111691583131</td><td>1·01094571509183</td></tr><tr><td>2</td><td>0·31698729810778</td><td>0·65036500052554</td><td>0·89220013824676</td></tr><tr><td>3</td><td>-0·18301270189222</td><td>-0·19093441556797</td><td>-0·03957502623564</td></tr><tr><td>4</td><td></td><td>-0·12083220831036</td><td>-0·26450716736904</td></tr><tr><td>5</td><td></td><td>0·04981749973164</td><td>0·04361630047418</td></tr><tr><td>6</td><td></td><td></td><td>0·04650360107098</td></tr><tr><td>7</td><td></td><td></td><td>-0·01498698933036</td></tr></tbody></table>
explicitly written as

$$ \phi(x) = a_0 \phi(2x) + a_1 \phi(2x-1) + \cdots + a_{N-1} \phi(2x-N+1) $$

Writing this relationship for all integer values $x = j$, it can be seen that

(a) for $j < 0$, $\phi(j)$ is a linear combination of $\phi(l < j)$. Thus if $\phi(x) \in L^2(\mathbf{R})$, i.e. square integrable, then there exists some integer $j_1 < 0$ for which $\phi(l < j_1) = 0$. This means that $\phi(j_1) = 0$. By induction, it follows that $\phi(j) = 0$ for all $j < 0$.

(b) for $j > N-1$, $\phi(j)$ is a linear combination of $\phi(l > j)$. Thus, if $\phi(x) \in L^2(\mathbf{R})$, then there exists some integer $j_2 > N-1$ for which $\phi(l > j_2) = 0$. This means that $\phi(j_2) = 0$. By induction, it follows that $\phi(j) = 0$ for all $j > N-1$.

Thus, the only remaining equations are

$$ \begin{aligned} \phi(0) &= a_0 \phi(0) \\ \phi(1) &= a_0 \phi(2) + a_1 \phi(1) + a_2 \phi(0) \\ \phi(2) &= a_0 \phi(4) + a_1 \phi(3) + a_2 \phi(2) + a_3 \phi(1) + a_4 \phi(0) \\ \vdots \\ \phi(N-2) &= a_{N-3} \phi(N-1) + a_{N-2} \phi(N-2) + a_{N-1} \phi(N-3) \\ \phi(N-1) &= a_{N-1} \phi(N-1) \end{aligned} $$

In matrix form this becomes

$$ \begin{bmatrix} a_0 & 0 & 0 & \cdots & 0 & 0 & 0 \\ a_2 & a_1 & a_0 & \cdots & 0 & 0 & 0 \\ a_4 & a_3 & a_2 & \cdots & 0 & 0 & 0 \\ \cdots & \cdots & \cdots & \cdots & \cdots & \cdots & \cdots \\ 0 & 0 & 0 & \cdots & a_{N-3} & a_{N-4} & a_{N-5} \\ 0 & 0 & 0 & \cdots & a_{N-1} & a_{N-2} & a_{N-3} \\ 0 & 0 & 0 & \cdots & 0 & 0 & a_{N-1} \end{bmatrix} \begin{bmatrix} \phi(0) \\ \phi(1) \\ \phi(2) \\ \cdots \\ \phi(N-3) \\ \phi(N-2) \\ \phi(N-1) \end{bmatrix} = \begin{bmatrix} \phi(0) \\ \phi(1) \\ \phi(2) \\ \cdots \\ \phi(N-3) \\ \phi(N-2) \\ \phi(N-1) \end{bmatrix} $$

or

$$ M\Phi = \Phi $$

Thus, the vector of integer values of the scaling function, $\Phi$, is the eigenvector of $M$ corresponding to the eigenvalue 1. As in all eigenvalue problems, the solution to the system

$$ (M - I)\Phi = 0 \qquad (27) $$

is not unique, and so a normalizing condition is required in order to determine a unique eigenvector. In this case the normalizing condition arises from equation (18). Some calculation shows that

$$ \sum_{k=-\infty}^{\infty} \phi(i) = 1, \quad i \in \mathbb{Z} \qquad (28) $$

Hence, the values of the scaling function at the integers are given by the solution to equation (27) normalized by equation (28).
Now that $\phi(x)$ is known at the integer values of $x$, the values of $\phi(x)$ at the half integers can be found from equation (16):

$$ \phi\left(\frac{x}{2}\right) = \sum_{k=-\infty}^{\infty} a_k \phi(x-k) \quad (29) $$

This process is repeated as many times as necessary to find the values of $\phi(x)$ at all dyadic points $\{i/2^n; i \in \mathbb{Z}\}$.

### 3.4. An example: the Daubechies 4 coefficient wavelet system

This section demonstrates how the results of the preceding discussion may be used to construct the Daubechies 4 coefficient (D4) scaling function and wavelet.

Consider first the construction of the filter coefficients. To solve equations (19), (21) and (26), it is easiest to include the redundant equation for $l=0$ in equation (26), which is linear, and exclude the non-linear equation for $l=1$ in equation (21). Thus,

$$ 
\begin{aligned}
a_0 + a_1 + a_2 + a_3 &= 2 \\
a_0^2 + a_1^2 + a_2^2 + a_3^2 &= 2 \\
a_0 - a_1 + a_2 - a_3 &= 0 \\
-a_1 + 2a_2 - 3a_3 &= 0
\end{aligned}
 $$

From the linear equations,

$$ a_0 = a_1 - \frac{1}{2}, \quad a_2 = \frac{1}{2} - a_1, \quad a_3 = 1 - a_1 $$

Substituting for $a_0, a_2$ and $a_3$ in the quadratic equation and solving gives

$$ a_1 = \frac{3 \pm \sqrt{3}}{4} $$

Thus,

$$ a_0 = \frac{1 \pm \sqrt{3}}{4}, \quad a_2 = \frac{3 \mp \sqrt{3}}{4}, \quad a_3 = \frac{1 \mp \sqrt{3}}{4} $$

Clearly, one set of solutions is the antithesis of the other. If one leads to the scaling function $\phi(x)$, then the other will lead to $\phi(-x)$. Both functions satisfy the conditions of normality, orthonormality and exact representation of a linear function. The solution that will be adopted is

$$ a_0 = \frac{1 + \sqrt{3}}{4}, \quad a_1 = \frac{3 + \sqrt{3}}{4}, \quad a_2 = \frac{3 - \sqrt{3}}{4}, \quad a_3 = \frac{1 - \sqrt{3}}{4} $$

The values of the scaling function at the integer points are then given by equation (27):

$$ \begin{bmatrix} a_0 - 1 & 0 & 0 & 0 \\ a_2 & a_1 - 1 & a_0 & 0 \\ 0 & a_3 & a_2 - 1 & a_1 \\ 0 & 0 & 0 & a_3 - 1 \end{bmatrix} \begin{bmatrix} \phi(0) \\ \phi(1) \\ \phi(2) \\ \phi(3) \end{bmatrix} = \begin{bmatrix} 0 \\ 0 \\ 0 \\ 0 \end{bmatrix} $$
from which

$$ \phi(0) = 0, \quad \phi(3) = 0, \quad \phi(2) = \frac{1 - \sqrt{3}}{1 + \sqrt{3}} \phi(1) $$

Substituting in equation (28) gives

$$ \phi(1) = \frac{1 + \sqrt{3}}{2} \quad \text{and} \quad \phi(2) = \frac{1 - \sqrt{3}}{2} $$

Equation (29) then gives the values of $\phi(x)$ at the half integers:

$$ \phi\left(\frac{1}{2}\right) = a_0 \phi(1) = \frac{2 + \sqrt{3}}{4} $$

$$ \phi\left(\frac{3}{2}\right) = a_1 \phi(2) + a_2 \phi(1) = 0 $$

$$ \phi\left(\frac{5}{2}\right) = a_3 \phi(2) + \frac{2 - \sqrt{3}}{4} $$

All other half integer values are zero, following an argument similar to that presented in Section 3.3 for the integer values of $\phi(x)$. Repeating this process for all dyadic points, it will be found that the scaling function is zero outside the interval $[0, 3]$, i.e. it has *compact support*.

## 4. CLASSIFICATION OF WAVELET BASES

There are many families of orthogonal wavelets that have been constructed in $L^2(\mathbf{R})$. Their differences lie in their localization in physical and frequency space, in their continuity and hence, in their differentiability. In 1982, Stromberg⁴ constructed wavelets which have $C^k$ continuity, where $k$ is arbitrary but finite. These wavelets have *exponential decay*,

$$ |\psi(x)| \le Ce^{-\gamma|x|} $$

In 1985, Meyer⁵ developed wavelets which, in *frequency space*, are both compactly supported and have $C^k$ continuity ($k$ is arbitrary and may be $\infty$). These wavelets are well localized in physical space too, showing decay faster than any inverse polynomial but not showing exponential decay. They are also $C^\infty$ continuous in physical space. In 1987, Tchamitchian constructed the first examples of *biorthogonal* wavelets. Soon after, Battle⁸ and Lemarié⁹ independently constructed identical families of orthonormal wavelets based on spline functions. These wavelets, like the Stromberg wavelets, have exponential decay (the decay rate decreases as $k$ increases) and are $C^k$ continuous. In 1986, Mallat⁶ and Meyer¹²,¹³ developed the theory of *multiresolution analysis*, which provides a general framework for the construction of wavelets.

The families of wavelets described so far exhibit inverse polynomial or exponential decay in physical space but they do not have compact support. Daubechies⁷ in 1988 provides the first general theory for the derivation of compactly supported wavelets. These wavelets are strictly localized in physical space while also retaining good localization in frequency space. Unlike the wavelets of Meyer, Battle and Lemarié, however, these wavelets do not exhibit symmetry.
## 5. MULTIRESOLUTION ANALYSIS USING THE MALLAT TRANSFORM

The Mallat transform¹¹ provides a simple means of transforming data from one level of resolution, *m*, to the next coarser level of resolution, *m* − 1. The inverse Mallat transform is a transformation from the coarser level, *m* − 1, back to the finer level, *m*.

### 5.1. Multiresolution decomposition

Multiresolution decomposition takes the expansion coefficients of the approximation, $P_m f$, to a function at scale *m* and decomposes them into

(1) the expansion coefficients, $c_{m-1,k}$, of the approximation, $P_{m-1} f$, at the next coarser scale, *m* − 1, and

(2) the expansion coefficients, $d_{m-1,k}$, of the detail component, $Q_{m-1} f = P_m f - P_{m-1} f$, at the next coarser scale.

Consider a function *f*. Let $P_m f$ denote the projection of *f* onto the subspace $V_m$ and $Q_m f$ denote the projection of *f* onto the subspace $W_m$. Thus,

$$P_m f = \sum_{k=-\infty}^{\infty} c_{m,k} \phi_{m,k}(x), \quad c_{m,k} = \langle f, \phi_{m,k} \rangle$$

and

$$Q_m f = \sum_{k=-\infty}^{\infty} d_{m,k} \psi_{m,k}(x), \quad d_{m,k} = \langle f, \psi_{m,k} \rangle$$

Since $W_{m-1}$ is the orthogonal complement of $V_{m-1}$ in $V_m$,

$$P_{m-1} f = P_m f - Q_{m-1} f$$

Substituting this in

$$c_{m-1,k} = \langle P_{m-1} f, \phi_{m-1,k} \rangle$$

leads to the following result:

$$c_{m-1,k} = \frac{1}{\sqrt{2}} \sum_{j=-\infty}^{\infty} c_{m,j} a_{j-2k} \quad (30)$$

Similarly, it can be shown that

$$d_{m-1,k} = \frac{1}{\sqrt{2}} \sum_{j=-\infty}^{\infty} c_{m,j} (-1)^j a_{N-1-j+2k} \quad (31)$$

Equations (30) and (31) form the basis of the Mallat transform algorithm described in Section 5.3.

### 5.2. Multiresolution reconstruction

Multiresolution reconstruction uses the scaling function and wavelet coefficients, $c_{m-1,k}$ and $d_{m-1,k}$, at level of resolution *m* − 1 to reconstruct the scaling function coefficients, $c_{m,k}$, at the next finer level of resolution, *m*.
Since $W_{m-1}$ is the orthogonal complement of $V_{m-1}$ in $V_m$,

$$P_m f = P_{m-1} f + Q_{m-1} f$$

Substituting this in

$$c_{m,k} = \langle P_m f, \phi_{m,k} \rangle$$

leads to

$$c_{m,k} = \frac{1}{\sqrt{2}} \sum_{j=-\infty}^{\infty} c_{m-1,j} a_{k-2j} + \frac{1}{\sqrt{2}} \sum_{j=-\infty}^{\infty} d_{m-1,j} (-1)^k a_{N-1-k+2j} \quad (32)$$

Equation (32) forms the basis of the inverse Mallat transform algorithm described in Section 5.3.

### 5.3. The Mallat transform and inverse transform algorithms

The results of Sections 5.1 and 5.2 have been expressed as a simple algorithm by Mallat.¹¹ The Mallat transform algorithm implements the decomposition process as follows.

Consider a string of data, $c_{m,k}$, of finite length, $n$, which represents the approximation, $P_m f$, to a function. For convenience, suppose that this data is periodic with period $n$. Then the matrix form of equation (30) is

$$ \begin{bmatrix} c_{m-1,0} \\ \times \\ c_{m-1,1} \\ \times \\ c_{m-1,2} \\ \times \\ \vdots \\ \vdots \\ c_{m-1,n/2-1} \\ \times \end{bmatrix} = \frac{1}{\sqrt{2}} \begin{bmatrix} a_0 & a_1 & a_2 & a_3 & \cdots & a_{N-1} & \cdots & 0 \\ 0 & a_0 & a_1 & a_2 & \cdots & a_{N-2} & \cdots & 0 \\ 0 & 0 & a_0 & a_1 & \cdots & a_{N-3} & \cdots & 0 \\ 0 & 0 & 0 & a_0 & \cdots & a_{N-4} & \cdots & 0 \\ \vdots & \vdots & \vdots & \vdots & \ddots & \vdots & \ddots & \vdots \\ 0 & 0 & 0 & 0 & \cdots & \cdots & a_{N-1} & \cdots \\ a_{N-1} & 0 & 0 & 0 & \cdots & \cdots & a_{N-2} & \cdots \\ \vdots & \vdots & \vdots & \vdots & \ddots & \ddots & \vdots & \ddots \\ a_2 & a_3 & a_4 & a_5 & \cdots & 0 & \cdots & a_1 \\ a_1 & a_2 & a_3 & a_4 & \cdots & 0 & \cdots & a_0 \end{bmatrix} \begin{bmatrix} c_{m,0} \\ c_{m,1} \\ c_{m,2} \\ c_{m,3} \\ \vdots \\ \vdots \\ c_{m,n-2} \\ c_{m,n-1} \end{bmatrix} $$

in which $\times$ represents information of no value. The effect of the periodicity is simply a wrap around of the filter coefficients at the bottom left corner of the matrix. Note that the matrix of filter coefficients is a circulant matrix, so that the right-hand side of the equation is effectively a convolution of the discrete filter

$$\tilde{h} = \frac{1}{\sqrt{2}} [a_0, 0, 0, 0, \dots, 0, a_{N-1}, \dots, a_2, a_1]^T$$
with the data $c_{m,k}$. Only every other sample of the result need be kept; this process is called decimation.

A similar process gives the expansion coefficients, $d_{m-1,k}$, of the detail which is lost in reducing the resolution of the data. The matrix form of equation (31) is

$$ \begin{bmatrix} d_{m-1,0} \\ \times \\ d_{m-1,1} \\ \times \\ d_{m-1,2} \\ \times \\ \vdots \\ \times \\ d_{m-1,n/2-1} \\ \times \end{bmatrix} = $$

$$ \frac{1}{\sqrt{2}} \begin{bmatrix} a_{N-1} & -a_{N-2} & a_{N-3} & -a_{N-4} & \cdots & -a_0 & \cdots & 0 \\ 0 & a_{N-1} & -a_{N-2} & a_{N-3} & \cdots & a_1 & \cdots & 0 \\ 0 & 0 & a_{N-1} & -a_{N-2} & \cdots & -a_2 & \cdots & 0 \\ 0 & 0 & 0 & a_{N-1} & \cdots & a_3 & \cdots & 0 \\ \vdots & \vdots & \vdots & \vdots & \ddots & \vdots & \ddots & \vdots \\ 0 & 0 & 0 & 0 & \cdots & 0 & \cdots & -a_0 \\ -a_0 & 0 & 0 & 0 & \cdots & 0 & \cdots & a_1 \\ \vdots & \vdots & \vdots & \vdots & \ddots & \vdots & \ddots & \vdots \\ a_{N-3} & -a_{N-4} & a_{N-5} & -a_{N-6} & \cdots & 0 & \cdots & -a_{N-2} \\ -a_{N-2} & a_{N-3} & -a_{N-4} & a_{N-5} & \cdots & 0 & \cdots & a_{N-1} \end{bmatrix} \begin{bmatrix} c_{m,0} \\ c_{m,1} \\ c_{m,2} \\ c_{m,3} \\ \vdots \\ c_{m,n-2} \\ c_{m,n-1} \end{bmatrix} $$

Once again, the matrix of filter coefficients is a circulant matrix, so that the right-hand side of the equation is effectively a convolution of the discrete filter

$$ \tilde{g} = \frac{1}{\sqrt{2}} [a_{N-1}, 0, 0, 0, \dots, 0, -a_0, \dots, a_{N-3}, -a_{N-2}]^T $$

with the data $c_{m,k}$.

Thus, the Mallat transform (decomposition) algorithm for the data, $c_{m,k}$, is

(a) convolve data with $\tilde{h}$ and decimate to get $c_{m-1,k}$,

(b) convolve data with $\tilde{g}$ and decimate to get $d_{m-1,k}$.

The inverse Mallat transform implements the reconstruction process. The matrix form of equation (32) is
$$ \begin{bmatrix} c_{m,0} \\ c_{m,1} \\ c_{m,2} \\ c_{m,3} \\ \vdots \\ \vdots \\ \vdots \\ c_{m,n-2} \\ c_{m,n-1} \end{bmatrix} = $$

$$ \frac{1}{\sqrt{2}} \begin{bmatrix} a_0 & 0 & 0 & 0 & \cdots & a_{N-1} & \cdots & a_1 \\ a_1 & a_0 & 0 & 0 & \cdots & 0 & \cdots & a_2 \\ a_2 & a_1 & a_0 & 0 & \cdots & 0 & \cdots & a_3 \\ a_3 & a_2 & a_1 & a_0 & \cdots & 0 & \cdots & a_4 \\ \vdots & \vdots & \vdots & \vdots & \ddots & \vdots & \ddots & \vdots \\ a_{N-2} & a_{N-3} & a_{N-4} & a_{N-5} & \cdots & \cdots & \cdots & a_{N-1} \\ a_{N-1} & a_{N-2} & a_{N-3} & a_{N-4} & \cdots & \cdots & \cdots & 0 \\ \vdots & \vdots & \vdots & \vdots & \ddots & \ddots & \ddots & \vdots \\ 0 & 0 & 0 & 0 & \cdots & a_{N-3} & \cdots & 0 \\ 0 & 0 & 0 & 0 & \cdots & a_{N-2} & \cdots & a_0 \end{bmatrix} \begin{bmatrix} c_{m-1,0} \\ 0 \\ c_{m-1,1} \\ 0 \\ c_{m-1,2} \\ 0 \\ \vdots \\ c_{m-1,n/2-1} \\ 0 \end{bmatrix} + $$

$$ \frac{1}{\sqrt{2}} \begin{bmatrix} a_{N-1} & 0 & 0 & 0 & \cdots & -a_0 & \cdots & -a_{N-2} \\ -a_{N-2} & a_{N-1} & 0 & 0 & \cdots & 0 & \cdots & a_{N-3} \\ a_{N-3} & -a_{N-2} & a_{N-1} & 0 & \cdots & 0 & \cdots & -a_{N-4} \\ -a_{N-4} & a_{N-3} & -a_{N-2} & a_{N-1} & \cdots & 0 & \cdots & a_{N-5} \\ \vdots & \vdots & \vdots & \vdots & \ddots & \vdots & \ddots & \vdots \\ a_1 & -a_2 & a_3 & -a_4 & \cdots & \cdots & \cdots & -a_0 \\ -a_0 & a_1 & -a_2 & a_3 & \cdots & \cdots & \cdots & 0 \\ \vdots & \vdots & \vdots & \vdots & \ddots & \vdots & \ddots & \vdots \\ 0 & 0 & 0 & 0 & \cdots & -a_2 & \cdots & 0 \\ 0 & 0 & 0 & 0 & \cdots & a_1 & \cdots & a_{N-1} \end{bmatrix} \begin{bmatrix} d_{m-1,0} \\ 0 \\ d_{m-1,1} \\ 0 \\ d_{m-1,2} \\ 0 \\ \vdots \\ d_{m-1,n/2-1} \\ 0 \end{bmatrix} $$

Thus, the inverse Mallat transform (reconstruction) algorithm is

(i) insert a zero between every sample in $c_{m-1,k}$ and $d_{m-1,k}$,

(ii) convolve $c_{m-1,k}$ with the filter $h$,

(iii) convolve $d_{m-1,k}$ with the filter $g$,

(iv) add, to get the original data, $c_{m,k}$,
where the filters *h* and *g* are given by

$$h = \frac{1}{\sqrt{2}} [a_0, a_1, a_2, a_3, \dots, a_{N-2}, a_{N-1}, \dots, 0, 0]^T$$

and

$$g = \frac{1}{\sqrt{2}} [a_{N-1}, -a_{N-2}, a_{N-3}, -a_{N-4}, \dots, a_1, -a_0, \dots, 0, 0]^T.$$

## 6. DATA ANALYSIS USING WAVELETS

One of the most promising applications of wavelets is in the analysis of engineering data. By nature of their construction, a scaling function is a low pass filter and a wavelet is a high pass filter. Figure 4 shows the frequency responses of the Battle-Lemarié scaling function, $\phi(x)$, and the Battle-Lemarié wavelet, $\psi(x)$, shown in Figure 3(g) and 3(h).

The effect of dilation by a factor of 2 in physical space is a corresponding scaling of the frequency axis. Thus, while $\psi(x)$ passes frequencies in the range [0.5, 1.0] Hz, $\psi(2x)$ passes frequencies in the range [1.0, 2.0] Hz. Furthermore, the good localization of the Battle-Lemarié wavelet in physical space means that the wavelet decomposition of a function will capture not only the frequencies contained in the function, but also the locations at which frequencies occur. Consider, for example, Figure 5(a), which shows a 24 Hz signal which has constant amplitude in the interval $x = [0, 0.5]$ and a 12 Hz modulated amplitude in the interval $x = [0.5, 1.0]$. Thus, the left half interval contains the frequency 24 Hz while the right half interval contains the frequencies 12 and 36 Hz. Figure 5(b)-(d) show the results of a wavelet decomposition at scales 4, 5 and 6, respectively. These clearly show the presence of a frequency in the left half interval at scale 5 ([16.0, 32.0] Hz) and the presence of two frequencies in the right half interval: one at scale 4 ([8.0, 16.0] Hz) and the other at scale 6 ([32.0, 64.0] Hz). This example clearly illustrates the time-frequency decomposition characteristics of wavelets. It also demonstrates that data may be compressed by discarding the wavelet components which are below a predetermined threshold.

### 6.1. Application to images

The most obvious two-dimensional application of wavelets is in image processing. An image may be regarded as a discrete sampling of a function of two variables, $f(x, y) \in [0, dx] \times [0, dy]$.

<figure><img src="image_5.png" /><figcaption>Figure 4. (a) Frequency response of the Battle-Lemarié scaling function. (b) Frequency response of the Battle-Lemarié wavelet</figcaption></figure>
<figure><img src="image_6.png" /><figcaption>Figure 5. (a) Original signal. (b) Wavelet coefficients of signal at scale 4. (c) Wavelet coefficients of signal at scale 5. (d) Wavelet coefficients of signal at scale 6</figcaption></figure>

Although by no means essential, it will be assumed here that $f(x, y)$ is periodic with periods $dx$ and $dy$ in $x$ and $y$. With this assumption, the Mallat transform algorithm of Section 5.3 may be conveniently applied to the decomposition of images.

The two-dimensional wavelet decomposition of the function, $f(x, y)$, is a tensor product:

$$f(x, y) = \sum_{j=-\infty}^{\infty} \sum_{k=-\infty}^{\infty} \sum_{l=-\infty}^{\infty} d_{j,k,l} \psi_{j,k}(x) \psi_{j,l}(y) \quad (33)$$

The first step in decomposing an image is to recognize that the image can be decomposed into its *scaling function* components at the finest scale, $m$, with $m$ being determined by the size of the image in pixels $(p, q)$:

$$f(x, y) \approx P_m f(x, y) = \sum_{k=-\infty}^{\infty} \sum_{l=-\infty}^{\infty} c_{m,k,l} \phi_{m,k}(x) \phi_{m,l}(y) \quad (34)$$

Since $m$ is the finest scale, this scaling function decomposition must result in $p \times q$ coefficients $c_{m,k,l}$, the maximum number permitted by the pixel data. It is assumed here that $p$ and $q$ are powers of 2, as implicitly demanded by the Mallat transform algorithm. From the assumption of periodicity and from equation (34), $m$ is given by

$$p = 2^m dx, \quad q = 2^m dy, \quad p, q \in \mathbb{Z}$$

The actual magnitudes of the assumed periods, $dx$ and $dy$, do not affect the resulting decomposition, but they must be consistent with the assumption that $p$ and $q$ are powers of 2 and their ratio must be consistent with the ratio $p/q$, e.g. for a 256 × 512 image, $dx$ might be 1 and $dy$ might be 2, so that $m = 8$.
In practice, the scaling function decomposition does not usually need to be carried out since, at the finest scale, the scaling function approaches the delta function, i.e.

$$ \phi_{m,k}(x) \rightarrow 2^{m/2} \delta(2^m x - k), \quad \text{as } m \rightarrow \infty $$

and this is the basis for a discrete sampling of $f(x, y)$. The projection of $f$ onto the space of delta functions is

$$ P_m f(x, y) = \sum_{k=-\infty}^{\infty} \sum_{l=-\infty}^{\infty} f_{k,l} 2^m \delta(2^m x - k) 2^m \delta(2^m y - l) \quad (35) $$

where

$f = f(2^{-mL}, 2^{-mL+1})$

<figure><img src="image_7.png" /><figcaption>Figure 6. Successive stages in the wavelet decomposition of an image, showing the averaged image and the detail components</figcaption></figure>
From equations (34) and (35), it is seen that for sufficiently large *m* (which is usually the case for images)

$$c_{m,k,l} \approx 2^m f_{k,l} = 2^m f(2^{-m}k, 2^{-m}l)$$

The $c_{m,k,l}$ may now be decomposed by repeatedly applying the two-dimensional Mallat transform algorithm, first to each of the columns of pixel values and then to each of the rows. After decomposition, it may be necessary to renormalize the pixel values if they lie outside the usual range [0, 255]. Figure 6 shows an image and the first three stages in the decomposition of the image. At each stage, the result is an averaged image and several detail components. For example, the first level of decomposition (Figure 6(b)) contains an average image corresponding to the basis

$$\phi_{m-1,k}(x) \phi_{m-1,l}(y)$$

and three detail components corresponding to the bases

$$\phi_{m-1,k}(x) \psi_{m-1,l}(y), \quad \psi_{m-1,k}(x) \phi_{m-1,l}(y) \quad \text{and} \quad \psi_{m-1,k}(x) \psi_{m-1,l}(y)$$

The average image may be further decomposed into four quadrants, while the detail components corresponding to the bases

$$\phi_{m-1,k}(x) \psi_{m-1,l}(y) \quad \text{and} \quad \psi_{m-1,k}(x) \phi_{m-1,l}(y)$$

may be decomposed into two halves. This gives the second level of decomposition (Figure 6(c)).

If this procedure is continued, the result would be a finite scale expansion similar to equation (33). On the other hand, it may be sufficient only to further decompose the average image at each stage, depending on the application in mind. As in the one-dimensional example, it is easy to see that efficient compression of the image may be obtained by storing only those components of the expansion which contain a significant amount of information and discarding the remaining components.

## 6.2. Application to stiffness matrices

The two-dimensional wavelet decomposition algorithm may be applied to any matrix. For example, a single stage wavelet decomposition produces an average matrix of $\frac{1}{4}$ the size of the original matrix, much like the average image in the top left-hand corner of Figure 6(b). The average matrix captures many of the properties of the original matrix. Consider, for example, the eigenvalues and eigenfrequencies of the stiffness matrix

$$ \begin{bmatrix} 1 & -1 & 0 & 0 & \dots & 0 & 0 \\ -1 & 2 & -1 & 0 & \dots & 0 & 0 \\ 0 & -1 & 2 & -1 & \dots & 0 & 0 \\ 0 & 0 & -1 & 2 & \dots & 0 & 0 \\ \vdots & \vdots & \vdots & \vdots & \ddots & \vdots & \vdots \\ 0 & 0 & 0 & 0 & \dots & 2 & -1 \\ 0 & 0 & 0 & 0 & \dots & -1 & 1 \end{bmatrix}_{64 \times 64} $$

Figure 7 compares the eigenvalues (and hence the eigenfrequencies) of the original matrix with those of the average matrices of $\frac{1}{4}$ and $\frac{1}{16}$ the size derived using the Haar, D4, D6 and D20 wavelet
<figure><img src="image_8.png" /><figcaption>Figure 7. Eigenvalues of a stiffness matrix and the wavelet reduced matrices using (a) Haar, (b) D4, (c) D6 and (d) D20</figcaption></figure>

systems. The figure shows that the eigenvalues of the average matrices are very similar to the lower eigenvalues of the original matrix. This result is to be expected since the average matrices are obtained by the low pass filtering action of the scaling function.

These results have significant implications in modal analysis, where engineers are often interested in determining the lowest eigenfrequencies since these are potentially the most damaging to a structure. The seismic damage to a structure can be catastrophic if its fundamental frequency is of the same order as the frequency of the earthquake. In this situation, therefore, it would be of interest to determine the lowest eigenvalue of the stiffness matrix and perhaps the eigenvalues corresponding to the second and third modes. The structure can have many more higher modes which are less important. Thus, by applying several stages of the wavelet decomposition to the stiffness matrix, good estimates of the most critical frequencies may be obtained with relatively little work.

## 7. SOLUTION OF PARTIAL DIFFERENTIAL EQUATIONS USING WAVELETS

The hierarchical nature of wavelets makes them a particularly appealing technique for the solution of partial differential equations. In the wavelet-Galerkin method,¹⁴⁻²⁰ the solution to a PDE is expressed as an expansion of scaling functions. Figure 8 shows how a function may be represented in terms of the Daubechies D6 scaling function at scale 1. The expansion coefficients represent the solution in wavelet space, i.e. they give the magnitude of each scaling function component of the solution.

The expansion coefficients may be calculated by solving a system of linear equations which contain the *connection coefficients*.¹⁵ The system of equations may be solved very efficiently using
<figure><img src="image_9.png" /><figcaption>Figure 8. Scaling function components of a given function u (solid line)</figcaption></figure>

the fast Fourier transform once it is recognized that it consists of one or more discrete convolutions. Full details of this procedure are given for one-dimensional PDEs in Reference 19 and for two-dimensional PDEs in Reference 20. Figure 9(a) shows a wavelet solution to Poisson's equation $u_{,xx} + u_{,yy} = f$ with Dirichlet boundary conditions. In order to facilitate comparison with the exact solution, the right-hand side $f$ was chosen to be that corresponding to

$$u = \sin\left(\frac{\pi}{6}x\right) \sin\left(\frac{7\pi}{4}x\right) \sin\left(\frac{3\pi}{4}y\right) \sin\left(\frac{5\pi}{4}y\right)$$

Figure 9(b) illustrates this exact solution.

Figure 10(a) shows the convergence of wavelet solutions to a one-dimensional problem as compared with the solution for a three-point finite difference scheme. Results were obtained using

<figure><img src="image_10.png" /><figcaption>Figure 9. (a) A wavelet solution to Poisson's equation (b) Exact solution</figcaption></figure>
<figure><img src="image_11.png" /><figcaption>Figure 10. (a) Convergence of Daubechies D6, D8, D10 and D12 wavelet solutions compared with finite difference solution (indicated by solid line) on a mesh with *n* unknowns (b) Accuracy vs. cost comparison of wavelet and finite difference solutions</figcaption></figure>

the Daubechies D6, D8, D10 and D12 scaling functions. The rate of convergence of the wavelet solutions increases dramatically as the order of the wavelet system increases, with almost no change in the computation time. The finite difference solution is slightly faster than the wavelet technique for a given mesh size since the latter involves the overhead of transforming to and from wavelet space. However, the overall accuracy vs. cost performance of the wavelet method is still far superior to the finite difference method (Figure 10(b)).

ACKNOWLEDGEMENT

This work was funded by a grant from NTT DATA Communications Systems Corporation, Kajima Corporation and Shimizu Corporation to the Intelligent Engineering Systems Laboratory, Massachusetts Institute of Technology.

REFERENCES

1. J. Morlet, G. Arens, I. Fourgeau and D. Giard, 'Wave propagation and sampling theory', *Geophys.*, **47**, 203-236 (1982).

2. J. Morlet, 'Sampling theory and wave propagation', in C. H. Chen (ed.), *NATO ASI Series, Issues in Acoustic Signal/Image Processing and Recognition*, Vol. I, Springer, Berlin, 1983, pp. 233-261.

3. A. Grossmann and J. Morlet, 'Decomposition of Hardy functions into square integrable wavelets of constant shape', *SIAM J. Math. Anal.*, **15**, 723-736 (1984).

4. J. O. Stromberg, 'A modified Franklin system and higher order spline systems on $R^n$ as unconditional bases for Hardy spaces', *Conf. in Honour of A. Zygmund*, Vol. II, *Wadsworth Math. Ser.*, 1982, pp. 475-493.

5. Y. Meyer, 'Principe d'incertitude, bases hilbertiennes et algèbres d'opérateurs', *Seminaire Bourbaki*, No. 662, 1985.

6. S. G. Mallat, 'Multiresolution approximation and wavelets', *Preprint, GRASP Lab., Department of Computer and Information Science, University of Pennsylvania*, 1986.

7. I. Daubechies, 'Orthonormal bases of compactly supported wavelets', *Commun. Pure Appl. Math.*, **41**, 909-996 (1988).

8. G. Battle, 'A block spin construction of ondelettes. Part I: Lemarie functions', *Comm. Math. Phys.*, **110**, 601-615 (1987).

9. P. G. Lemarie, 'Une nouvelle base d'ondettes de $L^2(R)$', *J. Math. Pures Appl.*, **67**, 227-236 (1988).

10. I. Daubechies, *Ten Lectures on Wavelets, CBMS-NSF Regional Conference Series, SIAM*, Philadelphia, 1992.

11. S. G. Mallat, 'A theory for multiresolution signal decomposition: the wavelet representation', *Commun. Pure Appl. Math.*, **41**(7), 674-693 (1988).

12. Y. Meyer, 'Ondettes et functions splines', *Lectures given at the University of Torino, Italy*, 1986.

13. Y. Meyer, 'Ondettes, fonctions splines et analyses graduees', *Seminaire EDP, Ecole Polytechnique, Paris, France*, 1986.
14. R. Glowinski, W. Lawton, M. Ravachol and E. Tenenbaum, 'Wavelet solution of linear and nonlinear elliptic, parabolic and hyperpolic problems in one space dimension', Proc. 9th Int. Conf. on Numerical Methods in Applied Sciences and Engineering, SIAM, Philadelphia, 1990.

15. A. Latto, H. Rasnikoff and E. Tenenbaum, 'The evaluation of connection coefficients of compactly supported wavelets', in Y. Maday (ed.), Proc. French-USA workshop on Wavelets and Turbulence, Princeton University, June 1991, Springer, Berlin, 1992.

16. A. Latto and E. Tenenbaum, 'Les ondelettes a support compact et la solution numerique de l'équation de Burgers', C. R. Acad. Sci. (Paris) 311, 903-909 (1990).

17. J. Weiss, 'Wavelets and the study of two dimensional turbulence', in Y. Maday (ed.), Proc. French-USA Workshop on Wavelets and Turbulence, Princeton University, June 1991, Springer, Berlin, 1992.

18. S. Qian and J. Weiss, 'Wavelets and the numerical solution of partial differential equations', J. Comp. Phys. 106(1), 155-175 (1993).

19. K. Amaratunga, J. R. Williams, S. Qian and J. Weiss, 'Wavelet-Galerkin solutions for one dimensional partial differential equations', Int. j. numer. methods eng., (1994) (to appear).

20. K. Amaratunga and J. R. Williams, 'Wavelet based Green's function approach to 2D PDEs', Eng. Comput., 10(4), 349-367 (1993).

21. W. H. Press, 'Wavelet transforms', Harvard-Smithsonian Center for Astrophysics Preprint No. 3184, 1991.

22. G. Beylkin, R. Coifman and V. Rokhlin, 'Fast wavelet transforms and numerical algorithms I', Comm. Pure Appl. Math., 44, 141-183 (1991).

23. C. K. Chui, *An Introduction to Wavelets*, Academic Press, New York, 1992.

24. C. K. Chui, *A Tutorial in Theory and Applications*, Academic Press, New York, 1992.

25. G. Strang, 'Wavelets and dilation equations: a brief introduction', SIAM Rev., 31(4), 614-627 (1989).

26. G. Strang, *Introduction to Applied Mathematics*, Wellesley-Cambridge Press, Wellesley, MA, 1986.