dwt.texi

This library is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY without ev

@cindex DWT, see wavelet transforms
@cindex wavelet transforms
@cindex transforms, wavelet

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 @file{gsl_wavelet.h} and @file{gsl_wavelet2d.h}.

@menu
* DWT Definitions::             
* DWT Initialization::          
* DWT Transform Functions::     
* DWT Examples::                
* DWT References::              
@end menu

@node DWT Definitions
@section Definitions
@cindex DWT, mathematical definition

The continuous wavelet transform and its inverse are defined by
the relations,
@iftex
@tex
\beforedisplay
$$
w(s, \tau) = \int_{-\infty}^\infty f(t) * \psi^*_{s,\tau}(t) dt
$$
\afterdisplay
@end tex
@end iftex
@ifinfo

@example
w(s,\tau) = \int f(t) * \psi^*_@{s,\tau@}(t) dt
@end example

@end ifinfo
@noindent
and,
@tex
\beforedisplay
$$
f(t) = \int_0^\infty ds \int_{-\infty}^\infty w(s, \tau) * \psi_{s,\tau}(t) d\tau
$$
\afterdisplay
@end tex
@ifinfo

@example
f(t) = \int \int_@{-\infty@}^\infty w(s, \tau) * \psi_@{s,\tau@}(t) d\tau ds
@end example

@end ifinfo
@noindent
where the basis functions @c{$\psi_{s,\tau}$}
@math{\psi_@{s,\tau@}} are obtained by scaling
and translation from a single function, referred to as the @dfn{mother
wavelet}.

The discrete version of the wavelet transform acts on equally-spaced
samples, with fixed scaling and translation steps (@math{s},
@math{\tau}).  The frequency and time axes are sampled @dfn{dyadically}
on scales of @math{2^j} through a level parameter @math{j}.
@c  The wavelet @math{\psi}
@c  can be expressed in terms of a scaling function @math{\varphi},
@c
@c  @tex
@c  \beforedisplay
@c  $$
@c  \psi(2^{j-1},t) = \sum_{k=0}^{2^j-1} g_j(k) * \bar{\varphi}(2^j t-k)
@c  $$
@c  \afterdisplay
@c  @end tex
@c  @ifinfo
@c  @example
@c  \psi(2^@{j-1@},t) = \sum_@{k=0@}^@{2^j-1@} g_j(k) * \bar@{\varphi@}(2^j t-k)
@c  @end example
@c  @end ifinfo
@c  @noindent
@c  and
@c
@c  @tex
@c  \beforedisplay
@c  $$
@c  \varphi(2^{j-1},t) = \sum_{k=0}^{2^j-1} h_j(k) * \bar{\varphi}(2^j t-k)
@c  $$
@c  \afterdisplay
@c  @end tex
@c  @ifinfo
@c  @example
@c  \varphi(2^@{j-1@},t) = \sum_@{k=0@}^@{2^j-1@} h_j(k) * \bar@{\varphi@}(2^j t-k)
@c  @end example
@c  @end ifinfo
@c  @noindent
@c  The functions @math{\psi} and @math{\varphi} are related through the
@c  coefficients
@c  @c{$g_{n} = (-1)^n h_{L-1-n}$}
@c  @math{g_@{n@} = (-1)^n h_@{L-1-n@}}
@c  for @c{$n=0 \dots L-1$}
@c  @math{n=0 ... L-1},
@c  where @math{L} is the total number of coefficients.  The two sets of
@c  coefficients @math{h_j} and @math{g_i} define the scaling function and
@c the wavelet.  
The resulting family of functions @c{$\{\psi_{j,n}\}$}
@math{@{\psi_@{j,n@}@}}
constitutes an orthonormal
basis for square-integrable signals.  

The discrete wavelet transform is an @math{O(N)} algorithm, and is also
referred to as the @dfn{fast wavelet transform}.

@node DWT Initialization
@section Initialization
@cindex DWT initialization

The @code{gsl_wavelet} structure contains the filter coefficients
defining the wavelet and any associated offset parameters.

@deftypefun {gsl_wavelet *} gsl_wavelet_alloc (const gsl_wavelet_type * @var{T}, size_t @var{k})
This function allocates and initializes a wavelet object of type
@var{T}.  The parameter @var{k} 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.
@end deftypefun

The following wavelet types are implemented:

@deffn {Wavelet} gsl_wavelet_daubechies
@deffnx {Wavelet} gsl_wavelet_daubechies_centered
@cindex Daubechies wavelets
@cindex maximal phase, Daubechies wavelets
The is the Daubechies wavelet family of maximum phase with @math{k/2}
vanishing moments.  The implemented wavelets are 
@c{$k=4, 6, \dots, 20$}
@math{k=4, 6, @dots{}, 20}, with @var{k} even.
@end deffn

@deffn {Wavelet} gsl_wavelet_haar
@deffnx {Wavelet} gsl_wavelet_haar_centered
@cindex Haar wavelets
This is the Haar wavelet.  The only valid choice of @math{k} for the
Haar wavelet is @math{k=2}.
@end deffn

@deffn {Wavelet} gsl_wavelet_bspline
@deffnx {Wavelet} gsl_wavelet_bspline_centered
@cindex biorthogonal wavelets
@cindex B-spline wavelets
This is the biorthogonal B-spline wavelet family of order @math{(i,j)}.  
The implemented values of @math{k = 100*i + j} are 103, 105, 202, 204,
206, 208, 301, 303, 305 307, 309.
@end deffn

@noindent
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.

@deftypefun {const char *} gsl_wavelet_name (const gsl_wavelet * @var{w})
This function returns a pointer to the name of the wavelet family for
@var{w}.
@end deftypefun

@c  @deftypefun {void} gsl_wavelet_print (const gsl_wavelet * @var{w})
@c  This function prints the filter coefficients (@code{**h1}, @code{**g1}, @code{**h2}, @code{**g2}) of the wavelet object @var{w}.
@c  @end deftypefun

@deftypefun {void} gsl_wavelet_free (gsl_wavelet * @var{w})
This function frees the wavelet object @var{w}.
@end deftypefun

The @code{gsl_wavelet_workspace} structure contains scratch space of the
same size as the input data and is used to hold intermediate results
during the transform.

@deftypefun {gsl_wavelet_workspace *} gsl_wavelet_workspace_alloc (size_t @var{n})
This function allocates a workspace for the discrete wavelet transform.
To perform a one-dimensional transform on @var{n} elements, a workspace
of size @var{n} must be provided.  For two-dimensional transforms of
@var{n}-by-@var{n} matrices it is sufficient to allocate a workspace of
size @var{n}, since the transform operates on individual rows and
columns.
@end deftypefun

@deftypefun {void} gsl_wavelet_workspace_free (gsl_wavelet_workspace * @var{work})
This function frees the allocated workspace @var{work}.
@end deftypefun

@node DWT Transform Functions
@section Transform Functions

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.

@menu
* DWT in one dimension::        
* DWT in two dimension::        
@end menu

@node DWT in one dimension
@subsection Wavelet transforms in one dimension
@cindex DWT, one dimensional

@deftypefun int gsl_wavelet_transform (const gsl_wavelet * @var{w}, double * @var{data}, size_t @var{stride}, size_t @var{n}, gsl_wavelet_direction @var{dir}, gsl_wavelet_workspace * @var{work})
@deftypefunx int gsl_wavelet_transform_forward (const gsl_wavelet * @var{w}, double * @var{data}, size_t @var{stride}, size_t @var{n}, gsl_wavelet_workspace * @var{work})
@deftypefunx int gsl_wavelet_transform_inverse (const gsl_wavelet * @var{w}, double * @var{data}, size_t @var{stride}, size_t @var{n}, gsl_wavelet_workspace * @var{work})

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

For the forward transform, the elements of the original array are 
replaced by the discrete wavelet
transform @c{$f_i \rightarrow w_{j,k}$}
@math{f_i -> w_@{j,k@}} 
in a packed triangular storage layout,