diff options
Diffstat (limited to 'gsl-1.9/doc/dwt.texi')
| -rw-r--r-- | gsl-1.9/doc/dwt.texi | 449 |
1 files changed, 449 insertions, 0 deletions
diff --git a/gsl-1.9/doc/dwt.texi b/gsl-1.9/doc/dwt.texi new file mode 100644 index 0000000..fa80b19 --- /dev/null +++ b/gsl-1.9/doc/dwt.texi @@ -0,0 +1,449 @@ +@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, +where @var{j} is the index of the level +@c{$j = 0 \dots J-1$} +@math{j = 0 ... J-1} +and +@var{k} is the index of the coefficient within each level, +@c{$k = 0 \dots 2^j - 1$} +@math{k = 0 ... (2^j)-1}. +The total number of levels is @math{J = \log_2(n)}. The output data +has the following form, +@tex +\beforedisplay +$$ +(s_{-1,0}, d_{0,0}, d_{1,0}, d_{1,1}, d_{2,0},\cdots, d_{j,k},\cdots, d_{J-1,2^{J-1} - 1}) +$$ +\afterdisplay +@end tex +@ifinfo + +@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@}) +@end example + +@end ifinfo +@noindent +where the first element is the smoothing coefficient @c{$s_{-1,0}$} +@math{s_@{-1,0@}}, followed by the detail coefficients @c{$d_{j,k}$} +@math{d_@{j,k@}} for each level +@math{j}. The backward transform inverts these coefficients to obtain +the original data. + +These functions return a status of @code{GSL_SUCCESS} upon successful +completion. @code{GSL_EINVAL} is returned if @var{n} is not an integer +power of 2 or if insufficient workspace is provided. +@end deftypefun + +@node DWT in two dimension +@subsection Wavelet transforms in two dimension +@cindex DWT, two dimensional + +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 complete discrete wavelet +transform on the rows of the matrix, followed by a separate complete +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 on the +rows and columns of the matrix for each level of the transform. The +first level of the transform is applied to the matrix rows, and then to +the matrix columns. 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 +@file{gsl_wavelet2d.h}. + +@deftypefun {int} gsl_wavelet2d_transform (const gsl_wavelet * @var{w}, double * @var{data}, size_t @var{tda}, size_t @var{size1}, size_t @var{size2}, gsl_wavelet_direction @var{dir}, gsl_wavelet_workspace * @var{work}) +@deftypefunx {int} gsl_wavelet2d_transform_forward (const gsl_wavelet * @var{w}, double * @var{data}, size_t @var{tda}, size_t @var{size1}, size_t @var{size2}, gsl_wavelet_workspace * @var{work}) +@deftypefunx {int} gsl_wavelet2d_transform_inverse (const gsl_wavelet * @var{w}, double * @var{data}, size_t @var{tda}, size_t @var{size1}, size_t @var{size2}, gsl_wavelet_workspace * @var{work}) + +These functions compute two-dimensional in-place forward and inverse +discrete wavelet transforms in standard and non-standard forms on the +array @var{data} stored in row-major form with dimensions @var{size1} +and @var{size2} and physical row length @var{tda}. The dimensions must +be equal (square matrix) and are 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 the appropriate size must be provided. On exit, +the appropriate elements of the array @var{data} are replaced by their +two-dimensional wavelet transform. + +The functions return a status of @code{GSL_SUCCESS} upon successful +completion. @code{GSL_EINVAL} is returned if @var{size1} and +@var{size2} are not equal and integer powers of 2, or if insufficient +workspace is provided. +@end deftypefun + +@deftypefun {int} gsl_wavelet2d_transform_matrix (const gsl_wavelet * @var{w}, gsl_matrix * @var{m}, gsl_wavelet_direction @var{dir}, gsl_wavelet_workspace * @var{work}) +@deftypefunx {int} gsl_wavelet2d_transform_matrix_forward (const gsl_wavelet * @var{w}, gsl_matrix * @var{m}, gsl_wavelet_workspace * @var{work}) +@deftypefunx {int} gsl_wavelet2d_transform_matrix_inverse (const gsl_wavelet * @var{w}, gsl_matrix * @var{m}, gsl_wavelet_workspace * @var{work}) +These functions compute the two-dimensional in-place wavelet transform +on a matrix @var{a}. +@end deftypefun + +@deftypefun {int} gsl_wavelet2d_nstransform (const gsl_wavelet * @var{w}, double * @var{data}, size_t @var{tda}, size_t @var{size1}, size_t @var{size2}, gsl_wavelet_direction @var{dir}, gsl_wavelet_workspace * @var{work}) +@deftypefunx {int} gsl_wavelet2d_nstransform_forward (const gsl_wavelet * @var{w}, double * @var{data}, size_t @var{tda}, size_t @var{size1}, size_t @var{size2}, gsl_wavelet_workspace * @var{work}) +@deftypefunx {int} gsl_wavelet2d_nstransform_inverse (const gsl_wavelet * @var{w}, double * @var{data}, size_t @var{tda}, size_t @var{size1}, size_t @var{size2}, gsl_wavelet_workspace * @var{work}) +These functions compute the two-dimensional wavelet transform in +non-standard form. +@end deftypefun + +@deftypefun {int} gsl_wavelet2d_nstransform_matrix (const gsl_wavelet * @var{w}, gsl_matrix * @var{m}, gsl_wavelet_direction @var{dir}, gsl_wavelet_workspace * @var{work}) +@deftypefunx {int} gsl_wavelet2d_nstransform_matrix_forward (const gsl_wavelet * @var{w}, gsl_matrix * @var{m}, gsl_wavelet_workspace * @var{work}) +@deftypefunx {int} gsl_wavelet2d_nstransform_matrix_inverse (const gsl_wavelet * @var{w}, gsl_matrix * @var{m}, gsl_wavelet_workspace * @var{work}) +These functions compute the non-standard form of the two-dimensional +in-place wavelet transform on a matrix @var{a}. +@end deftypefun + +@node DWT Examples +@section Examples + +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. + +@example +@verbatiminclude examples/dwt.c +@end example + +@noindent +The output can be used with the @sc{gnu} plotutils @code{graph} program, + +@example +$ ./a.out ecg.dat > dwt.dat +$ graph -T ps -x 0 256 32 -h 0.3 -a dwt.dat > dwt.ps +@end example + +@iftex +The graphs below show an original and compressed version of a sample ECG +recording from the MIT-BIH Arrhythmia Database, part of the PhysioNet +archive of public-domain of medical datasets. + +@sp 1 +@center @image{dwt-orig,3.4in} +@center @image{dwt-samp,3.4in} +@quotation +Original (upper) and wavelet-compressed (lower) ECG signals, using the +20 largest components of the Daubechies(4) discrete wavelet transform. +@end quotation +@end iftex + +@node DWT References +@section References and Further Reading + +The mathematical background to wavelet transforms is covered in the +original lectures by Daubechies, + +@itemize @asis +@item +Ingrid Daubechies. +Ten Lectures on Wavelets. +@cite{CBMS-NSF Regional Conference Series in Applied Mathematics} (1992), +SIAM, ISBN 0898712742. +@end itemize + +@noindent +An easy to read introduction to the subject with an emphasis on the +application of the wavelet transform in various branches of science is, + +@itemize @asis +@item +Paul S. Addison. @cite{The Illustrated Wavelet Transform Handbook}. +Institute of Physics Publishing (2002), ISBN 0750306920. +@end itemize + +@noindent +For extensive coverage of signal analysis by wavelets, wavelet packets +and local cosine bases see, + +@itemize @asis +@item +S. G. Mallat. @cite{A wavelet tour of signal processing} (Second +edition). Academic Press (1999), ISBN 012466606X. +@end itemize + +@noindent +The concept of multiresolution analysis underlying the wavelet transform +is described in, + +@itemize @asis +@item +S. G. Mallat. +Multiresolution Approximations and Wavelet Orthonormal Bases of L@math{^2}(R). +@cite{Transactions of the American Mathematical Society}, 315(1), 1989, 69--87. +@end itemize + +@itemize @asis +@item +S. G. Mallat. +A Theory for Multiresolution Signal Decomposition---The Wavelet Representation. +@cite{IEEE Transactions on Pattern Analysis and Machine Intelligence}, 11, 1989, +674--693. +@end itemize + +@noindent +The coefficients for the individual wavelet families implemented by the +library can be found in the following papers, + +@itemize @asis +@item +I. Daubechies. +Orthonormal Bases of Compactly Supported Wavelets. +@cite{Communications on Pure and Applied Mathematics}, 41 (1988) 909--996. +@end itemize + +@itemize @asis +@item +A. Cohen, I. Daubechies, and J.-C. Feauveau. +Biorthogonal Bases of Compactly Supported Wavelets. +@cite{Communications on Pure and Applied Mathematics}, 45 (1992) +485--560. +@end itemize + +@noindent +The PhysioNet archive of physiological datasets can be found online at +@uref{http://www.physionet.org/} and is described in the following +paper, + +@itemize @asis +@item +Goldberger et al. +PhysioBank, PhysioToolkit, and PhysioNet: Components +of a New Research Resource for Complex Physiologic +Signals. +@cite{Circulation} 101(23):e215-e220 2000. +@end itemize |