Go to the first, previous, next, last section, table of contents.
This chapter describes functions for performing Discrete Wavelet Transforms (DWTs). The library includes wavelets for real data in both one and two dimensions. The wavelet functions are declared in the header files 'gsl_wavelet.h' and 'gsl_wavelet2d.h'.
The continuous wavelet transform and its inverse are defined by the relations,
w(s,\tau) = \int f(t) * \psi^*_{s,\tau}(t) dt
and,
f(t) = \int \int_{-\infty}^\infty w(s, \tau) * \psi_{s,\tau}(t) d\tau ds
where the basis functions \psi_{s,\tau} are obtained by scaling and translation from a single function, referred to as the mother wavelet.
The discrete version of the wavelet transform acts on evenly sampled data, with fixed scaling and translation steps (s, \tau). The frequency and time axes are sampled dyadically on scales of 2^j through a level parameter j. The resulting family of functions {\psi_{j,n}} constitutes an orthonormal basis for square-integrable signals.
The discrete wavelet transform is an O(N) algorithm, and is also referred to as the fast wavelet transform.
The gsl_wavelet structure contains the filter coefficients
defining the wavelet and associated offset parameters (for wavelets with
centered support).
The following wavelet types are implemented:
The centered forms of the wavelets align the coefficients of the various sub-bands on edges. Thus the resulting visualization of the coefficients of the wavelet transform in the phase plane is easier to understand.
The gsl_wavelet_workspace structure contains scratch space of the
same size as the input data, for holding intermediate results during the
transform.
This sections describes the actual functions performing the discrete wavelet transform. Note that the transforms use periodic boundary conditions. If the signal is not periodic in the sample length then spurious coefficients will appear at the beginning and end of each level of the transform.
These functions compute in-place forward and inverse discrete wavelet
transforms of length n with stride stride on the array
data. The length of the transform n is restricted to powers
of two. For the transform version of the function the argument
dir can be either forward (+1) or backward
(-1). A workspace work of length n must be provided.
For the forward transform, the elements of the original array are replaced by the discrete wavelet transform f_i -> w_{j,k} in a packed triangular storage layout, where j is the index of the level j = 0 ... J-1 and k is the index of the coefficient within each level, k = 0 ... (2^j)-1. The total number of levels is J = \log_2(n). The output data has the following form,
(s_{-1,0}, d_{0,0}, d_{1,0}, d_{1,1}, d_{2,0}, ...,
d_{j,k}, ..., d_{J-1,2^{J-1}-1})
where the first element is the smoothing coefficient s_{-1,0}, followed by the detail coefficients d_{j,k} for each level j. The backward transform inverts these coefficients to obtain the original data.
These functions return a status of GSL_SUCCESS upon successful
completion. GSL_EINVAL is returned if n is not an integer
power of 2 or if insufficient workspace is provided.
The library provides functions to perform two-dimensional discrete wavelet transforms on square matrices. The matrix dimensions must be an integer power of two. There are two possible orderings of the rows and columns in the two-dimensional wavelet transform, referred to as the "standard" and "non-standard" forms.
The "standard" transform performs a discrete wavelet transform on all rows of the matrix, followed by a separate discrete wavelet transform on the columns of the resulting row-transformed matrix. This procedure uses the same ordering as a two-dimensional fourier transform.
The "non-standard" transform is performed in interleaved passes of the each level of the transform on the rows and columns of the matrix. The first level of the transform is carried out on the matrix rows, and then the columns of the partially row-transformed data. This procedure is then repeated across the rows and columns of the data for the subsequent levels of the transform, until the full discrete wavelet transform is complete. The non-standard form of the discrete wavelet transform is typically used in image analysis.
The functions described in this section are declared in the header file 'gsl_wavelet2d.h'.
These functions compute two-dimensional in-place forward and inverse
discrete wavelet transforms in standard and non-standard forms on the
array data stored in row-major form with dimensions size1
and size2 and physical row length tda. The dimensions must
be equal (square matrix) and are restricted to powers of two. For the
transform version of the function the argument dir can be
either forward (+1) or backward (-1). A
workspace work of the appropriate size must be provided. On exit,
the appropriate elements of the array data are replaced by their
two-dimensional wavelet transform.
The functions return a status of GSL_SUCCESS upon successful
completion. GSL_EINVAL is returned if size1 and
size2 are not equal and integer powers of 2, or if insufficient
workspace is provided.
The following program demonstrates the use of the one-dimensional wavelet transform functions. It computes an approximation to an input signal (of length 256) using the 20 largest components of the wavelet transform, while setting the others to zero.
#include <stdio.h>
#include <math.h>
#include <gsl/gsl_sort.h>
#include <gsl/gsl_wavelet.h>
int
main (int argc, char **argv)
{
int i, n = 256, nc = 20;
double *data = malloc (n * sizeof (double));
double *abscoeff = malloc (n * sizeof (double));
size_t *p = malloc (n * sizeof (size_t));
FILE *f = fopen (argv[1], "r");
for (i = 0; i < n; i++)
{
fscanf (f, "%lg", &data[i]);
}
fclose (f);
{
gsl_wavelet *w = gsl_wavelet_alloc (gsl_wavelet_daubechies, 4);
gsl_wavelet_workspace *work = gsl_wavelet_workspace_alloc (n);
gsl_wavelet_transform_forward (w, data, 1, n, work);
for (i = 0; i < n; i++)
{
abscoeff[i] = fabs (data[i]);
}
gsl_sort_index (p, abscoeff, 1, n);
for (i = 0; (i + nc) < n; i++)
data[p[i]] = 0;
gsl_wavelet_transform_inverse (w, data, 1, n, work);
}
for (i = 0; i < n; i++)
{
printf ("%g\n", data[i]);
}
}
The output can be used with the GNU plotutils graph program,
$ ./a.out ecg.dat > dwt.dat $ graph -T ps -x 0 256 32 -h 0.3 -a dwt.dat > dwt.ps
The mathematical background to wavelet transforms is covered in the original lectures by Daubechies,
An easy to read introduction to the subject with an emphasis on the application of the wavelet transform in various branches of science is,
For extensive coverage of signal analysis by wavelets, wavelet packets and local cosine bases see,
The concept of multiresolution analysis underlying the wavelet transform is described in,
The coefficients for the individual wavelet families implemented by the library can be found in the following papers,
The PhysioNet archive of physiological datasets can be found online at http://www.physionet.org/ and is described in the following paper,
Go to the first, previous, next, last section, table of contents.