File: gsl-ref_30.html

package info (click to toggle)
gsl-ref-html 1.6-1
links: PTS
area: main
in suites: sarge
size: 1,504 kB
ctags: 3,558
sloc: makefile: 36
file content (568 lines) | stat: -rw-r--r-- 18,366 bytes
<HTML>
<HEAD>
<!-- This HTML file has been created by texi2html 1.54+ (gsl)
     from ../gsl-ref.texi -->

<TITLE>GNU Scientific Library -- Reference Manual - Wavelet Transforms</TITLE>
<!-- <LINK rel="stylesheet" title="Default Style Sheet" href="/css/texinfo.css" type="text/css"> -->
<link href="gsl-ref_31.html" rel=Next>
<link href="gsl-ref_29.html" rel=Previous>
<link href="gsl-ref_toc.html" rel=ToC>

</HEAD>
<BODY>
<p>Go to the <A HREF="gsl-ref_1.html">first</A>, <A HREF="gsl-ref_29.html">previous</A>, <A HREF="gsl-ref_31.html">next</A>, <A HREF="gsl-ref_50.html">last</A> section, <A HREF="gsl-ref_toc.html">table of contents</A>.
<P><HR><P>


<H1><A NAME="SEC416" HREF="gsl-ref_toc.html#TOC416">Wavelet Transforms</A></H1>
<P>
<A NAME="IDX2013"></A>
<A NAME="IDX2014"></A>
<A NAME="IDX2015"></A>
<A NAME="IDX2016"></A>
<A NAME="IDX2017"></A>

</P>
<P>
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 <TT>'gsl_wavelet.h'</TT> and <TT>'gsl_wavelet2d.h'</TT>.

</P>



<H2><A NAME="SEC417" HREF="gsl-ref_toc.html#TOC417">Definitions</A></H2>
<P>
<A NAME="IDX2018"></A>

</P>
<P>
The continuous wavelet transform and its inverse are defined by
the relations,

</P>

<PRE class="example">
w(s,\tau) = \int f(t) * \psi^*_{s,\tau}(t) dt
</PRE>

<P>
and,

</P>

<PRE class="example">
f(t) = \int \int_{-\infty}^\infty w(s, \tau) * \psi_{s,\tau}(t) d\tau ds
</PRE>

<P>
where the basis functions 
\psi_{s,\tau} are obtained by scaling
and translation from a single function, referred to as the <I>mother
wavelet</I>.

</P>
<P>
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 <I>dyadically</I> 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.  

</P>
<P>
The discrete wavelet transform is an O(N) algorithm, and is also
referred to as the <I>fast wavelet transform</I>.

</P>


<H2><A NAME="SEC418" HREF="gsl-ref_toc.html#TOC418">Initialization</A></H2>
<P>
<A NAME="IDX2019"></A>

</P>
<P>
The <CODE>gsl_wavelet</CODE> structure contains the filter coefficients
defining the wavelet and associated offset parameters (for wavelets with
centered support).  

</P>
<P>
<DL>
<DT><U>Function:</U> gsl_wavelet * <B>gsl_wavelet_alloc</B> <I>(const gsl_wavelet_type * <VAR>T</VAR>, size_t <VAR>k</VAR>)</I>
<DD><A NAME="IDX2020"></A>
This function allocates and initializes a wavelet object of type
<VAR>T</VAR>.  The parameter <VAR>k</VAR> selects the specific member of the
wavelet family.  A null pointer is returned if insufficient memory is
available or if a unsupported member is selected.
</DL>

</P>
<P>
The following wavelet types are implemented:

</P>
<P>
<DL>
<DT><U>Wavelet:</U> <B>gsl_wavelet_daubechies</B>
<DD><A NAME="IDX2021"></A>
<DT><U>Wavelet:</U> <B>gsl_wavelet_daubechies_centered</B>
<DD><A NAME="IDX2022"></A>
<A NAME="IDX2023"></A>
<A NAME="IDX2024"></A>
The is the Daubechies wavelet family of maximum phase with k/2
vanishing moments.  The implemented wavelets are 
k=4, 6, ..., 20, with <VAR>k</VAR> even.
</DL>

</P>
<P>
<DL>
<DT><U>Wavelet:</U> <B>gsl_wavelet_haar</B>
<DD><A NAME="IDX2025"></A>
<DT><U>Wavelet:</U> <B>gsl_wavelet_haar_centered</B>
<DD><A NAME="IDX2026"></A>
<A NAME="IDX2027"></A>
This is the Haar wavelet.  The only valid choice of k for the
Haar wavelet is k=2.
</DL>

</P>
<P>
<DL>
<DT><U>Wavelet:</U> <B>gsl_wavelet_bspline</B>
<DD><A NAME="IDX2028"></A>
<DT><U>Wavelet:</U> <B>gsl_wavelet_bspline_centered</B>
<DD><A NAME="IDX2029"></A>
<A NAME="IDX2030"></A>
<A NAME="IDX2031"></A>
This is the biorthogonal B-spline wavelet family of order (i,j).  
The implemented values of k = 100*i + j are 103, 105, 202, 204,
206, 208, 301, 303, 305 307, 309.
</DL>

</P>
<P>
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.

</P>
<P>
<DL>
<DT><U>Function:</U> const char * <B>gsl_wavelet_name</B> <I>(const gsl_wavelet * <VAR>w</VAR>)</I>
<DD><A NAME="IDX2032"></A>
This function returns a pointer to the name of the wavelet family for
<VAR>w</VAR>.
</DL>

</P>

<P>
<DL>
<DT><U>Function:</U> void <B>gsl_wavelet_free</B> <I>(gsl_wavelet * <VAR>w</VAR>)</I>
<DD><A NAME="IDX2033"></A>
This function frees the wavelet object <VAR>w</VAR>.
</DL>

</P>
<P>
The <CODE>gsl_wavelet_workspace</CODE> structure contains scratch space of the
same size as the input data, for holding intermediate results during the
transform.

</P>
<P>
<DL>
<DT><U>Function:</U> gsl_wavelet_workspace * <B>gsl_wavelet_workspace_alloc</B> <I>(size_t <VAR>n</VAR>)</I>
<DD><A NAME="IDX2034"></A>
This function allocates a workspace for the discrete wavelet transform.
To perform a one-dimensional transform on <VAR>n</VAR> elements, a workspace
of size <VAR>n</VAR> must be provided.  For two-dimensional transforms of
<VAR>n</VAR>-by-<VAR>n</VAR> matrices it is sufficient to allocate a workspace of
size <VAR>n</VAR>, since the transform operates on individual rows and
columns.
</DL>

</P>
<P>
<DL>
<DT><U>Function:</U> void <B>gsl_wavelet_workspace_free</B> <I>(gsl_wavelet_workspace * <VAR>workspace</VAR>)</I>
<DD><A NAME="IDX2035"></A>
This function frees the allocated workspace <VAR>workspace</VAR>.
</DL>

</P>


<H2><A NAME="SEC419" HREF="gsl-ref_toc.html#TOC419">Transform Functions</A></H2>

<P>
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.

</P>



<H3><A NAME="SEC420" HREF="gsl-ref_toc.html#TOC420">Wavelet transforms in one dimension</A></H3>
<P>
<A NAME="IDX2036"></A>

</P>
<P>
<DL>
<DT><U>Function:</U> int <B>gsl_dwt_transform</B> <I>(const gsl_wavelet * <VAR>w</VAR>, double * <VAR>data</VAR>, size_t <VAR>stride</VAR>, size_t <VAR>n</VAR>, gsl_wavelet_direction <VAR>dir</VAR>, gsl_wavelet_workspace * <VAR>work</VAR>)</I>
<DD><A NAME="IDX2037"></A>
<DT><U>Function:</U> int <B>gsl_dwt_transform_forward</B> <I>(const gsl_wavelet * <VAR>w</VAR>, double * <VAR>data</VAR>, size_t <VAR>stride</VAR>, size_t <VAR>n</VAR>, gsl_wavelet_workspace * <VAR>work</VAR>)</I>
<DD><A NAME="IDX2038"></A>
<DT><U>Function:</U> int <B>gsl_dwt_transform_inverse</B> <I>(const gsl_wavelet * <VAR>w</VAR>, double * <VAR>data</VAR>, size_t <VAR>stride</VAR>, size_t <VAR>n</VAR>, gsl_wavelet_workspace * <VAR>work</VAR>)</I>
<DD><A NAME="IDX2039"></A>

</P>
<P>
These functions compute in-place forward and inverse discrete wavelet
transforms of length <VAR>n</VAR> with stride <VAR>stride</VAR> on the array
<VAR>data</VAR>. The length of the transform <VAR>n</VAR> is restricted to powers
of two.  For the <CODE>transform</CODE> version of the function the argument
<VAR>dir</VAR> can be either <CODE>forward</CODE> (+1) or <CODE>backward</CODE>
(-1).  A workspace <VAR>work</VAR> of length <VAR>n</VAR> must be provided.

</P>
<P>
For the forward transform, the elements of the original array are 
replaced by the discrete wavelet
transform 
f_i -&#62; w_{j,k} 
in a packed triangular storage layout, 
where <VAR>j</VAR> is the index of the level 
j = 0 ... J-1
and
<VAR>k</VAR> 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,

</P>

<PRE class="example">
(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}) 
</PRE>

<P>
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.

</P>
<P>
These functions return a status of <CODE>GSL_SUCCESS</CODE> upon successful
completion.  <CODE>GSL_EINVAL</CODE> is returned if <VAR>n</VAR> is not an integer
power of 2 or if insufficient workspace is provided.
</DL>

</P>


<H3><A NAME="SEC421" HREF="gsl-ref_toc.html#TOC421">Wavelet transforms in two dimension</A></H3>
<P>
<A NAME="IDX2040"></A>

</P>
<P>
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.

</P>
<P>
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.

</P>
<P>
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.

</P>
<P>
The functions described in this section are declared in the header file
<TT>'gsl_wavelet2d.h'</TT>.

</P>
<P>
<DL>
<DT><U>Function:</U> int <B>gsl_wavelet2d_transform</B> <I>(const gsl_wavelet * <VAR>w</VAR>, double * <VAR>data</VAR>, size_t <VAR>tda</VAR>, size_t <VAR>size1</VAR>, size_t <VAR>size2</VAR>, gsl_wavelet_direction <VAR>dir</VAR>, gsl_wavelet_workspace * <VAR>work</VAR>)</I>
<DD><A NAME="IDX2041"></A>
<DT><U>Function:</U> int <B>gsl_wavelet2d_transform_forward</B> <I>(const gsl_wavelet * <VAR>w</VAR>, double * <VAR>data</VAR>, size_t <VAR>tda</VAR>, size_t <VAR>size1</VAR>, size_t <VAR>size2</VAR>, gsl_wavelet_workspace * <VAR>work</VAR>)</I>
<DD><A NAME="IDX2042"></A>
<DT><U>Function:</U> int <B>gsl_wavelet2d_transform_inverse</B> <I>(const gsl_wavelet * <VAR>w</VAR>, double * <VAR>data</VAR>, size_t <VAR>tda</VAR>, size_t <VAR>size1</VAR>, size_t <VAR>size2</VAR>, gsl_wavelet_workspace * <VAR>work</VAR>)</I>
<DD><A NAME="IDX2043"></A>

</P>
<P>
These functions compute two-dimensional in-place forward and inverse
discrete wavelet transforms in standard and non-standard forms on the
array <VAR>data</VAR> stored in row-major form with dimensions <VAR>size1</VAR>
and <VAR>size2</VAR> and physical row length <VAR>tda</VAR>.  The dimensions must
be equal (square matrix) and are restricted to powers of two.  For the
<CODE>transform</CODE> version of the function the argument <VAR>dir</VAR> can be
either <CODE>forward</CODE> (+1) or <CODE>backward</CODE> (-1).  A
workspace <VAR>work</VAR> of the appropriate size must be provided.  On exit,
the appropriate elements of the array <VAR>data</VAR> are replaced by their
two-dimensional wavelet transform.

</P>
<P>
The functions return a status of <CODE>GSL_SUCCESS</CODE> upon successful
completion.  <CODE>GSL_EINVAL</CODE> is returned if <VAR>size1</VAR> and
<VAR>size2</VAR> are not equal and integer powers of 2, or if insufficient
workspace is provided.
</DL>

</P>
<P>
<DL>
<DT><U>Function:</U> int <B>gsl_wavelet2d_transform_matrix</B> <I>(const gsl_wavelet * <VAR>w</VAR>, gsl_matrix * <VAR>m</VAR>, gsl_wavelet_direction <VAR>dir</VAR>, gsl_wavelet_workspace * <VAR>work</VAR>)</I>
<DD><A NAME="IDX2044"></A>
<DT><U>Function:</U> int <B>gsl_wavelet2d_transform_matrix_forward</B> <I>(const gsl_wavelet * <VAR>w</VAR>, gsl_matrix * <VAR>m</VAR>, gsl_wavelet_workspace * <VAR>work</VAR>)</I>
<DD><A NAME="IDX2045"></A>
<DT><U>Function:</U> int <B>gsl_wavelet2d_transform_matrix_inverse</B> <I>(const gsl_wavelet * <VAR>w</VAR>, gsl_matrix * <VAR>m</VAR>, gsl_wavelet_workspace * <VAR>work</VAR>)</I>
<DD><A NAME="IDX2046"></A>
These functions compute the two-dimensional in-place wavelet transform
on a matrix <VAR>a</VAR>.
</DL>

</P>
<P>
<DL>
<DT><U>Function:</U> int <B>gsl_wavelet2d_nstransform</B> <I>(const gsl_wavelet * <VAR>w</VAR>, double * <VAR>data</VAR>, size_t <VAR>tda</VAR>, size_t <VAR>size1</VAR>, size_t <VAR>size2</VAR>, gsl_wavelet_direction <VAR>dir</VAR>, gsl_wavelet_workspace * <VAR>work</VAR>)</I>
<DD><A NAME="IDX2047"></A>
<DT><U>Function:</U> int <B>gsl_wavelet2d_nstransform_forward</B> <I>(const gsl_wavelet * <VAR>w</VAR>, double * <VAR>data</VAR>, size_t <VAR>tda</VAR>, size_t <VAR>size1</VAR>, size_t <VAR>size2</VAR>, gsl_wavelet_workspace * <VAR>work</VAR>)</I>
<DD><A NAME="IDX2048"></A>
<DT><U>Function:</U> int <B>gsl_wavelet2d_nstransform_inverse</B> <I>(const gsl_wavelet * <VAR>w</VAR>, double * <VAR>data</VAR>, size_t <VAR>tda</VAR>, size_t <VAR>size1</VAR>, size_t <VAR>size2</VAR>, gsl_wavelet_workspace * <VAR>work</VAR>)</I>
<DD><A NAME="IDX2049"></A>
These functions compute the two-dimensional wavelet transform in
non-standard form.
</DL>

</P>
<P>
<DL>
<DT><U>Function:</U> int <B>gsl_wavelet2d_nstransform_matrix</B> <I>(const gsl_wavelet * <VAR>w</VAR>, gsl_matrix * <VAR>m</VAR>, gsl_wavelet_direction <VAR>dir</VAR>, gsl_wavelet_workspace * <VAR>work</VAR>)</I>
<DD><A NAME="IDX2050"></A>
<DT><U>Function:</U> int <B>gsl_wavelet2d_nstransform_matrix_forward</B> <I>(const gsl_wavelet * <VAR>w</VAR>, gsl_matrix * <VAR>m</VAR>, gsl_wavelet_workspace * <VAR>work</VAR>)</I>
<DD><A NAME="IDX2051"></A>
<DT><U>Function:</U> int <B>gsl_wavelet2d_nstransform_matrix_inverse</B> <I>(const gsl_wavelet * <VAR>w</VAR>, gsl_matrix * <VAR>m</VAR>, gsl_wavelet_workspace * <VAR>work</VAR>)</I>
<DD><A NAME="IDX2052"></A>
These functions compute the non-standard form of the two-dimensional
in-place wavelet transform on a matrix <VAR>a</VAR>.
</DL>

</P>


<H2><A NAME="SEC422" HREF="gsl-ref_toc.html#TOC422">Example</A></H2>

<P>
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.

</P>

<PRE class="example">
#include &#60;stdio.h&#62;
#include &#60;math.h&#62;
#include &#60;gsl/gsl_sort.h&#62;
#include &#60;gsl/gsl_wavelet.h&#62;

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 &#60; n; i++)
    {
      fscanf (f, "%lg", &#38;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 &#60; n; i++)
      {
        abscoeff[i] = fabs (data[i]);
      }

    gsl_sort_index (p, abscoeff, 1, n);

    for (i = 0; (i + nc) &#60; n; i++)
      data[p[i]] = 0;

    gsl_wavelet_transform_inverse (w, data, 1, n, work);
  }

  for (i = 0; i &#60; n; i++)
    {
      printf ("%g\n", data[i]);
    }
}
</PRE>

<P>
The output can be used with the GNU plotutils <CODE>graph</CODE> program,

</P>

<PRE class="example">
$ ./a.out ecg.dat &#62; dwt.dat
$ graph -T ps -x 0 256 32 -h 0.3 -a dwt.dat &#62; dwt.ps
</PRE>



<H2><A NAME="SEC423" HREF="gsl-ref_toc.html#TOC423">References and Further Reading</A></H2>

<P>
The mathematical background to wavelet transforms is covered in the
original lectures by Daubechies,

</P>

<UL class="itemize">
<LI>

Ingrid Daubechies.
Ten Lectures on Wavelets.
<CITE>CBMS-NSF Regional Conference Series in Applied Mathematics</CITE> (1992), 
SIAM, ISBN 0898712742.
</UL>

<P>
An easy to read introduction to the subject with an emphasis on the
application of the wavelet transform in various branches of science is,

</P>

<UL class="itemize">
<LI>

Paul S. Addison. <CITE>The Illustrated Wavelet Transform Handbook</CITE>.
Institute of Physics Publishing (2002), ISBN 0750306920.
</UL>

<P>
For extensive coverage of signal analysis by wavelets, wavelet packets
and local cosine bases see,

</P>

<UL class="itemize">
<LI>

St&eacute;phane Mallat.  <CITE>A wavelet tour of signal processing</CITE> (Second
edition). Academic Press (1999), ISBN 012466606X.
</UL>

<P>
The concept of multiresolution analysis underlying the wavelet transform
is described in,

</P>

<UL class="itemize">
<LI>

S. G. Mallat.
Multiresolution Approximations and Wavelet Orthonormal Bases of L^2(R).
<CITE>Transactions of the American Mathematical Society</CITE>, 315(1), 1989, 69--87.
</UL>


<UL class="itemize">
<LI>

S. G. Mallat.
A Theory for Multiresolution Signal Decomposition -- The Wavelet Representation.
<CITE>IEEE Transactions on Pattern Analysis and Machine Intelligence</CITE>, 11, 1989,
674--693. 
</UL>

<P>
The coefficients for the individual wavelet families implemented by the
library can be found in the following papers,

</P>

<UL class="itemize">
<LI>

I. Daubechies.
Orthonormal Bases of Compactly Supported Wavelets.
<CITE>Communications on Pure and Applied Mathematics</CITE>, 41 (1988) 909--996.
</UL>


<UL class="itemize">
<LI>

A. Cohen, I. Daubechies, and J.-C. Feauveau.
Biorthogonal Bases of Compactly Supported Wavelets.
<CITE>Communications on Pure and Applied Mathematics</CITE>, 45 (1992)
485--560.
</UL>

<P>
The PhysioNet archive of physiological datasets can be found online at
<A HREF="http://www.physionet.org/">http://www.physionet.org/</A> and is described in the following
paper,

</P>

<UL class="itemize">
<LI>

Goldberger et al.  
PhysioBank, PhysioToolkit, and PhysioNet: Components
of a New Research Resource for Complex Physiologic
Signals. 
<CITE>Circulation</CITE> 101(23):e215-e220 2000.
</UL>

<P><HR><P>
<p>Go to the <A HREF="gsl-ref_1.html">first</A>, <A HREF="gsl-ref_29.html">previous</A>, <A HREF="gsl-ref_31.html">next</A>, <A HREF="gsl-ref_50.html">last</A> section, <A HREF="gsl-ref_toc.html">table of contents</A>.
</BODY>
</HTML>