htswanalysis/MACS/lib/gsl/gsl-1.11/doc/dwt.texi

   1 @cindex DWT, see wavelet transforms
   2 @cindex wavelet transforms
   3 @cindex transforms, wavelet
   4
   5 This chapter describes functions for performing Discrete Wavelet
   6 Transforms (DWTs).  The library includes wavelets for real data in both
   7 one and two dimensions.  The wavelet functions are declared in the header
   8 files @file{gsl_wavelet.h} and @file{gsl_wavelet2d.h}.
   9
  10 @menu
  11 * DWT Definitions::
  12 * DWT Initialization::
  13 * DWT Transform Functions::
  14 * DWT Examples::
  15 * DWT References::
  16 @end menu
  17
  18 @node DWT Definitions
  19 @section Definitions
  20 @cindex DWT, mathematical definition
  21
  22 The continuous wavelet transform and its inverse are defined by
  23 the relations,
  24 @iftex
  25 @tex
  26 \beforedisplay
  27 $$
  28 w(s, \tau) = \int_{-\infty}^\infty f(t) * \psi^*_{s,\tau}(t) dt
  29 $$
  30 \afterdisplay
  31 @end tex
  32 @end iftex
  33 @ifinfo
  34
  35 @example
  36 w(s,\tau) = \int f(t) * \psi^*_@{s,\tau@}(t) dt
  37 @end example
  38
  39 @end ifinfo
  40 @noindent
  41 and,
  42 @tex
  43 \beforedisplay
  44 $$
  45 f(t) = \int_0^\infty ds \int_{-\infty}^\infty w(s, \tau) * \psi_{s,\tau}(t) d\tau
  46 $$
  47 \afterdisplay
  48 @end tex
  49 @ifinfo
  50
  51 @example
  52 f(t) = \int \int_@{-\infty@}^\infty w(s, \tau) * \psi_@{s,\tau@}(t) d\tau ds
  53 @end example
  54
  55 @end ifinfo
  56 @noindent
  57 where the basis functions @c{$\psi_{s,\tau}$}
  58 @math{\psi_@{s,\tau@}} are obtained by scaling
  59 and translation from a single function, referred to as the @dfn{mother
  60 wavelet}.
  61
  62 The discrete version of the wavelet transform acts on equally-spaced
  63 samples, with fixed scaling and translation steps (@math{s},
  64 @math{\tau}).  The frequency and time axes are sampled @dfn{dyadically}
  65 on scales of @math{2^j} through a level parameter @math{j}.
  66 @c  The wavelet @math{\psi}
  67 @c  can be expressed in terms of a scaling function @math{\varphi},
  68 @c
  69 @c  @tex
  70 @c  \beforedisplay
  71 @c  $$
  72 @c  \psi(2^{j-1},t) = \sum_{k=0}^{2^j-1} g_j(k) * \bar{\varphi}(2^j t-k)
  73 @c  $$
  74 @c  \afterdisplay
  75 @c  @end tex
  76 @c  @ifinfo
  77 @c  @example
  78 @c  \psi(2^@{j-1@},t) = \sum_@{k=0@}^@{2^j-1@} g_j(k) * \bar@{\varphi@}(2^j t-k)
  79 @c  @end example
  80 @c  @end ifinfo
  81 @c  @noindent
  82 @c  and
  83 @c
  84 @c  @tex
  85 @c  \beforedisplay
  86 @c  $$
  87 @c  \varphi(2^{j-1},t) = \sum_{k=0}^{2^j-1} h_j(k) * \bar{\varphi}(2^j t-k)
  88 @c  $$
  89 @c  \afterdisplay
  90 @c  @end tex
  91 @c  @ifinfo
  92 @c  @example
  93 @c  \varphi(2^@{j-1@},t) = \sum_@{k=0@}^@{2^j-1@} h_j(k) * \bar@{\varphi@}(2^j t-k)
  94 @c  @end example
  95 @c  @end ifinfo
  96 @c  @noindent
  97 @c  The functions @math{\psi} and @math{\varphi} are related through the
  98 @c  coefficients
  99 @c  @c{$g_{n} = (-1)^n h_{L-1-n}$}
 100 @c  @math{g_@{n@} = (-1)^n h_@{L-1-n@}}
 101 @c  for @c{$n=0 \dots L-1$}
 102 @c  @math{n=0 ... L-1},
 103 @c  where @math{L} is the total number of coefficients.  The two sets of
 104 @c  coefficients @math{h_j} and @math{g_i} define the scaling function and
 105 @c the wavelet.
 106 The resulting family of functions @c{$\{\psi_{j,n}\}$}
 107 @math{@{\psi_@{j,n@}@}}
 108 constitutes an orthonormal
 109 basis for square-integrable signals.
 110
 111 The discrete wavelet transform is an @math{O(N)} algorithm, and is also
 112 referred to as the @dfn{fast wavelet transform}.
 113
 114 @node DWT Initialization
 115 @section Initialization
 116 @cindex DWT initialization
 117
 118 The @code{gsl_wavelet} structure contains the filter coefficients
 119 defining the wavelet and any associated offset parameters.
 120
 121 @deftypefun {gsl_wavelet *} gsl_wavelet_alloc (const gsl_wavelet_type * @var{T}, size_t @var{k})
 122 This function allocates and initializes a wavelet object of type
 123 @var{T}.  The parameter @var{k} selects the specific member of the
 124 wavelet family.  A null pointer is returned if insufficient memory is
 125 available or if a unsupported member is selected.
 126 @end deftypefun
 127
 128 The following wavelet types are implemented:
 129
 130 @deffn {Wavelet} gsl_wavelet_daubechies
 131 @deffnx {Wavelet} gsl_wavelet_daubechies_centered
 132 @cindex Daubechies wavelets
 133 @cindex maximal phase, Daubechies wavelets
 134 The is the Daubechies wavelet family of maximum phase with @math{k/2}
 135 vanishing moments.  The implemented wavelets are
 136 @c{$k=4, 6, \dots, 20$}
 137 @math{k=4, 6, @dots{}, 20}, with @var{k} even.
 138 @end deffn
 139
 140 @deffn {Wavelet} gsl_wavelet_haar
 141 @deffnx {Wavelet} gsl_wavelet_haar_centered
 142 @cindex Haar wavelets
 143 This is the Haar wavelet.  The only valid choice of @math{k} for the
 144 Haar wavelet is @math{k=2}.
 145 @end deffn
 146
 147 @deffn {Wavelet} gsl_wavelet_bspline
 148 @deffnx {Wavelet} gsl_wavelet_bspline_centered
 149 @cindex biorthogonal wavelets
 150 @cindex B-spline wavelets
 151 This is the biorthogonal B-spline wavelet family of order @math{(i,j)}.
 152 The implemented values of @math{k = 100*i + j} are 103, 105, 202, 204,
 153 206, 208, 301, 303, 305 307, 309.
 154 @end deffn
 155
 156 @noindent
 157 The centered forms of the wavelets align the coefficients of the various
 158 sub-bands on edges.  Thus the resulting visualization of the
 159 coefficients of the wavelet transform in the phase plane is easier to
 160 understand.
 161
 162 @deftypefun {const char *} gsl_wavelet_name (const gsl_wavelet * @var{w})
 163 This function returns a pointer to the name of the wavelet family for
 164 @var{w}.
 165 @end deftypefun
 166
 167 @c  @deftypefun {void} gsl_wavelet_print (const gsl_wavelet * @var{w})
 168 @c  This function prints the filter coefficients (@code{**h1}, @code{**g1}, @code{**h2}, @code{**g2}) of the wavelet object @var{w}.
 169 @c  @end deftypefun
 170
 171 @deftypefun {void} gsl_wavelet_free (gsl_wavelet * @var{w})
 172 This function frees the wavelet object @var{w}.
 173 @end deftypefun
 174
 175 The @code{gsl_wavelet_workspace} structure contains scratch space of the
 176 same size as the input data and is used to hold intermediate results
 177 during the transform.
 178
 179 @deftypefun {gsl_wavelet_workspace *} gsl_wavelet_workspace_alloc (size_t @var{n})
 180 This function allocates a workspace for the discrete wavelet transform.
 181 To perform a one-dimensional transform on @var{n} elements, a workspace
 182 of size @var{n} must be provided.  For two-dimensional transforms of
 183 @var{n}-by-@var{n} matrices it is sufficient to allocate a workspace of
 184 size @var{n}, since the transform operates on individual rows and
 185 columns.
 186 @end deftypefun
 187
 188 @deftypefun {void} gsl_wavelet_workspace_free (gsl_wavelet_workspace * @var{work})
 189 This function frees the allocated workspace @var{work}.
 190 @end deftypefun
 191
 192 @node DWT Transform Functions
 193 @section Transform Functions
 194
 195 This sections describes the actual functions performing the discrete
 196 wavelet transform.  Note that the transforms use periodic boundary
 197 conditions.  If the signal is not periodic in the sample length then
 198 spurious coefficients will appear at the beginning and end of each level
 199 of the transform.
 200
 201 @menu
 202 * DWT in one dimension::
 203 * DWT in two dimension::
 204 @end menu
 205
 206 @node DWT in one dimension
 207 @subsection Wavelet transforms in one dimension
 208 @cindex DWT, one dimensional
 209
 210 @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})
 211 @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})
 212 @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})
 213
 214 These functions compute in-place forward and inverse discrete wavelet
 215 transforms of length @var{n} with stride @var{stride} on the array
 216 @var{data}. The length of the transform @var{n} is restricted to powers
 217 of two.  For the @code{transform} version of the function the argument
 218 @var{dir} can be either @code{forward} (@math{+1}) or @code{backward}
 219 (@math{-1}).  A workspace @var{work} of length @var{n} must be provided.
 220
 221 For the forward transform, the elements of the original array are
 222 replaced by the discrete wavelet
 223 transform @c{$f_i \rightarrow w_{j,k}$}
 224 @math{f_i -> w_@{j,k@}}
 225 in a packed triangular storage layout,
 226 where @var{j} is the index of the level
 227 @c{$j = 0 \dots J-1$}
 228 @math{j = 0 ... J-1}
 229 and
 230 @var{k} is the index of the coefficient within each level,
 231 @c{$k = 0 \dots 2^j - 1$}
 232 @math{k = 0 ... (2^j)-1}.
 233 The total number of levels is @math{J = \log_2(n)}.  The output data
 234 has the following form,
 235 @tex
 236 \beforedisplay
 237 $$
 238 (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})
 239 $$
 240 \afterdisplay
 241 @end tex
 242 @ifinfo
 243
 244 @example
 245 (s_@{-1,0@}, d_@{0,0@}, d_@{1,0@}, d_@{1,1@}, d_@{2,0@}, ...,
 246   d_@{j,k@}, ..., d_@{J-1,2^@{J-1@}-1@})
 247 @end example
 248
 249 @end ifinfo
 250 @noindent
 251 where the first element is the smoothing coefficient @c{$s_{-1,0}$}
 252 @math{s_@{-1,0@}}, followed by the detail coefficients @c{$d_{j,k}$}
 253 @math{d_@{j,k@}} for each level
 254 @math{j}.  The backward transform inverts these coefficients to obtain
 255 the original data.
 256
 257 These functions return a status of @code{GSL_SUCCESS} upon successful
 258 completion.  @code{GSL_EINVAL} is returned if @var{n} is not an integer
 259 power of 2 or if insufficient workspace is provided.
 260 @end deftypefun
 261
 262 @node DWT in two dimension
 263 @subsection Wavelet transforms in two dimension
 264 @cindex DWT, two dimensional
 265
 266 The library provides functions to perform two-dimensional discrete
 267 wavelet transforms on square matrices.  The matrix dimensions must be an
 268 integer power of two.  There are two possible orderings of the rows and
 269 columns in the two-dimensional wavelet transform, referred to as the
 270 ``standard'' and ``non-standard'' forms.
 271
 272 The ``standard'' transform performs a complete discrete wavelet
 273 transform on the rows of the matrix, followed by a separate complete
 274 discrete wavelet transform on the columns of the resulting
 275 row-transformed matrix.  This procedure uses the same ordering as a
 276 two-dimensional fourier transform.
 277
 278 The ``non-standard'' transform is performed in interleaved passes on the
 279 rows and columns of the matrix for each level of the transform.  The
 280 first level of the transform is applied to the matrix rows, and then to
 281 the matrix columns.  This procedure is then repeated across the rows and
 282 columns of the data for the subsequent levels of the transform, until
 283 the full discrete wavelet transform is complete.  The non-standard form
 284 of the discrete wavelet transform is typically used in image analysis.
 285
 286 The functions described in this section are declared in the header file
 287 @file{gsl_wavelet2d.h}.
 288
 289 @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})
 290 @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})
 291 @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})
 292
 293 These functions compute two-dimensional in-place forward and inverse
 294 discrete wavelet transforms in standard and non-standard forms on the
 295 array @var{data} stored in row-major form with dimensions @var{size1}
 296 and @var{size2} and physical row length @var{tda}.  The dimensions must
 297 be equal (square matrix) and are restricted to powers of two.  For the
 298 @code{transform} version of the function the argument @var{dir} can be
 299 either @code{forward} (@math{+1}) or @code{backward} (@math{-1}).  A
 300 workspace @var{work} of the appropriate size must be provided.  On exit,
 301 the appropriate elements of the array @var{data} are replaced by their
 302 two-dimensional wavelet transform.
 303
 304 The functions return a status of @code{GSL_SUCCESS} upon successful
 305 completion.  @code{GSL_EINVAL} is returned if @var{size1} and
 306 @var{size2} are not equal and integer powers of 2, or if insufficient
 307 workspace is provided.
 308 @end deftypefun
 309
 310 @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})
 311 @deftypefunx {int} gsl_wavelet2d_transform_matrix_forward (const gsl_wavelet * @var{w}, gsl_matrix * @var{m}, gsl_wavelet_workspace * @var{work})
 312 @deftypefunx {int} gsl_wavelet2d_transform_matrix_inverse (const gsl_wavelet * @var{w}, gsl_matrix * @var{m}, gsl_wavelet_workspace * @var{work})
 313 These functions compute the two-dimensional in-place wavelet transform
 314 on a matrix @var{a}.
 315 @end deftypefun
 316
 317 @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})
 318 @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})
 319 @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})
 320 These functions compute the two-dimensional wavelet transform in
 321 non-standard form.
 322 @end deftypefun
 323
 324 @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})
 325 @deftypefunx {int} gsl_wavelet2d_nstransform_matrix_forward (const gsl_wavelet * @var{w}, gsl_matrix * @var{m}, gsl_wavelet_workspace * @var{work})
 326 @deftypefunx {int} gsl_wavelet2d_nstransform_matrix_inverse (const gsl_wavelet * @var{w}, gsl_matrix * @var{m}, gsl_wavelet_workspace * @var{work})
 327 These functions compute the non-standard form of the two-dimensional
 328 in-place wavelet transform on a matrix @var{a}.
 329 @end deftypefun
 330
 331 @node DWT Examples
 332 @section Examples
 333
 334 The following program demonstrates the use of the one-dimensional
 335 wavelet transform functions.  It computes an approximation to an input
 336 signal (of length 256) using the 20 largest components of the wavelet
 337 transform, while setting the others to zero.
 338
 339 @example
 340 @verbatiminclude examples/dwt.c
 341 @end example
 342
 343 @noindent
 344 The output can be used with the @sc{gnu} plotutils @code{graph} program,
 345
 346 @example
 347 $ ./a.out ecg.dat > dwt.dat
 348 $ graph -T ps -x 0 256 32 -h 0.3 -a dwt.dat > dwt.ps
 349 @end example
 350
 351 @iftex
 352 The graphs below show an original and compressed version of a sample ECG
 353 recording from the MIT-BIH Arrhythmia Database, part of the PhysioNet
 354 archive of public-domain of medical datasets.
 355
 356 @sp 1
 357 @center @image{dwt-orig,3.4in}
 358 @center @image{dwt-samp,3.4in}
 359 @quotation
 360 Original (upper) and wavelet-compressed (lower) ECG signals, using the
 361 20 largest components of the Daubechies(4) discrete wavelet transform.
 362 @end quotation
 363 @end iftex
 364
 365 @node DWT References
 366 @section References and Further Reading
 367
 368 The mathematical background to wavelet transforms is covered in the
 369 original lectures by Daubechies,
 370
 371 @itemize @asis
 372 @item
 373 Ingrid Daubechies.
 374 Ten Lectures on Wavelets.
 375 @cite{CBMS-NSF Regional Conference Series in Applied Mathematics} (1992),
 376 SIAM, ISBN 0898712742.
 377 @end itemize
 378
 379 @noindent
 380 An easy to read introduction to the subject with an emphasis on the
 381 application of the wavelet transform in various branches of science is,
 382
 383 @itemize @asis
 384 @item
 385 Paul S. Addison. @cite{The Illustrated Wavelet Transform Handbook}.
 386 Institute of Physics Publishing (2002), ISBN 0750306920.
 387 @end itemize
 388
 389 @noindent
 390 For extensive coverage of signal analysis by wavelets, wavelet packets
 391 and local cosine bases see,
 392
 393 @itemize @asis
 394 @item
 395 S. G. Mallat.  @cite{A wavelet tour of signal processing} (Second
 396 edition). Academic Press (1999), ISBN 012466606X.
 397 @end itemize
 398
 399 @noindent
 400 The concept of multiresolution analysis underlying the wavelet transform
 401 is described in,
 402
 403 @itemize @asis
 404 @item
 405 S. G. Mallat.
 406 Multiresolution Approximations and Wavelet Orthonormal Bases of L@math{^2}(R).
 407 @cite{Transactions of the American Mathematical Society}, 315(1), 1989, 69--87.
 408 @end itemize
 409
 410 @itemize @asis
 411 @item
 412 S. G. Mallat.
 413 A Theory for Multiresolution Signal Decomposition---The Wavelet Representation.
 414 @cite{IEEE Transactions on Pattern Analysis and Machine Intelligence}, 11, 1989,
 415 674--693.
 416 @end itemize
 417
 418 @noindent
 419 The coefficients for the individual wavelet families implemented by the
 420 library can be found in the following papers,
 421
 422 @itemize @asis
 423 @item
 424 I. Daubechies.
 425 Orthonormal Bases of Compactly Supported Wavelets.
 426 @cite{Communications on Pure and Applied Mathematics}, 41 (1988) 909--996.
 427 @end itemize
 428
 429 @itemize @asis
 430 @item
 431 A. Cohen, I. Daubechies, and J.-C. Feauveau.
 432 Biorthogonal Bases of Compactly Supported Wavelets.
 433 @cite{Communications on Pure and Applied Mathematics}, 45 (1992)
 434 485--560.
 435 @end itemize
 436
 437 @noindent
 438 The PhysioNet archive of physiological datasets can be found online at
 439 @uref{http://www.physionet.org/} and is described in the following
 440 paper,
 441
 442 @itemize @asis
 443 @item
 444 Goldberger et al.
 445 PhysioBank, PhysioToolkit, and PhysioNet: Components
 446 of a New Research Resource for Complex Physiologic
 447 Signals.
 448 @cite{Circulation} 101(23):e215-e220 2000.
 449 @end itemize