Main Page   Namespace List   Class Hierarchy   Compound List   File List   Compound Members   File Members  

Daubechies Class Template Reference

Daubechies D4 wavelet transform (D4 denotes four coefficients). More...

#include <daub.h>

Inheritance diagram for Daubechies::

liftbase List of all members.

Public Methods

void forwardStep (T &a, const int n)
 Forward Daubechies D4 transform step. More...

void forwardStepRev (T &a, const int n)
 Forward Daubechies D4 transform step, where the locations for the high and low pass filters are reversed. More...

void inverseStep (T &a, const int n)
 Inverse Daubechies D4 transform. More...

 Daubechies ()
 Initialize the filter constants used in the Daubechies D4 transform. More...

void forwardTrans (T &ts, int N)
 Forward Daubechies D4 transform. More...

void inverseTrans (T &coef, int N)
 Inverse Daubechies D4 transform. More...


Protected Methods

void predict (T &vec, int N, transDirection direction)
 Predict step, to be defined by the subclass. More...

virtual void update (T &vec, int N, transDirection direction)
 Update step, to be defined by the subclass. More...


Private Attributes

double h0
 forward transform scaling coefficients. More...

double h1
 forward transform scaling coefficients. More...

double h2
 forward transform scaling coefficients. More...

double h3
 forward transform scaling coefficients. More...

double g0
 forward transform wave coefficients. More...

double g1
 forward transform wave coefficients. More...

double g2
 forward transform wave coefficients. More...

double g3
 forward transform wave coefficients. More...

double Ih0
double Ih1
double Ih2
double Ih3
double Ig0
double Ig1
double Ig2
double Ig3

Detailed Description

template<class T> class Daubechies

Daubechies D4 wavelet transform (D4 denotes four coefficients).

I have to confess up front that the comment here does not even come close to describing wavelet algorithms and the Daubechies D4 algorithm in particular. I don't think that it can be described in anything less than a journal article or perhaps a book. I even have to apologize for the notation I use to describe the algorithm, which is barely adequate. But explaining the correct notation would take a fair amount of space as well. This comment really represents some notes that I wrote up as I implemented the code. If you are unfamiliar with wavelets I suggest that you look at the bearcave.com web pages and at the wavelet literature. I have yet to see a really good reference on wavelets for the software developer. The best book I can recommend is Ripples in Mathematics by Jensen and Cour-Harbo.

All wavelet algorithms have two components, a wavelet function and a scaling function. These are sometime also referred to as high pass and low pass filters respectively.

The wavelet function is passed two or more samples and calculates a wavelet coefficient. In the case of the Haar wavelet this is

  coefi = oddi - eveni
  or 
  coefi = 0.5 * (oddi - eveni)
  

depending on the version of the Haar algorithm used.

The scaling function produces a smoother version of the original data. In the case of the Haar wavelet algorithm this is an average of two adjacent elements.

The Daubechies D4 wavelet algorithm also has a wavelet and a scaling function. The coefficients for the scaling function are denoted as hi and the wavelet coefficients are gi.

Mathematicians like to talk about wavelets in terms of a wavelet algorithm applied to an infinite data set. In this case one step of the forward transform can be expressed as the infinite matrix of wavelet coefficients represented below multiplied by the infinite signal vector.

     ai = ...h0,h1,h2,h3, 0, 0, 0, 0, 0, 0, 0, ...   si
     ci = ...g0,g1,g2,g3, 0, 0, 0, 0, 0, 0, 0, ...   si+1
   ai+1 = ...0, 0, h0,h1,h2,h3, 0, 0, 0, 0, 0, ...   si+2
   ci+1 = ...0, 0, g0,g1,g2,g3, 0, 0, 0, 0, 0, ...   si+3
   ai+2 = ...0, 0, 0, 0, h0,h1,h2,h3, 0, 0, 0, ...   si+4
   ci+2 = ...0, 0, 0, 0, g0,g1,g2,g3, 0, 0, 0, ...   si+5
   ai+3 = ...0, 0, 0, 0, 0, 0, h0,h1,h2,h3, 0, ...   si+6
   ci+3 = ...0, 0, 0, 0, 0, 0, g0,g1,g2,g3, 0, ...   si+7
  

The dot product (inner product) of the infinite vector and a row of the matrix produces either a smoother version of the signal (ai) or a wavelet coefficient (ci).

In an ordered wavelet transform, the smoothed (ai) are stored in the first half of an n element array region. The wavelet coefficients (ci) are stored in the second half the n element region. The algorithm is recursive. The smoothed values become the input to the next step.

The transpose of the forward transform matrix above is used to calculate an inverse transform step. Here the dot product is formed from the result of the forward transform and an inverse transform matrix row.

      si = ...h2,g2,h0,g0, 0, 0, 0, 0, 0, 0, 0, ...  ai
    si+1 = ...h3,g3,h1,g1, 0, 0, 0, 0, 0, 0, 0, ...  ci
    si+2 = ...0, 0, h2,g2,h0,g0, 0, 0, 0, 0, 0, ...  ai+1
    si+3 = ...0, 0, h3,g3,h1,g1, 0, 0, 0, 0, 0, ...  ci+1
    si+4 = ...0, 0, 0, 0, h2,g2,h0,g0, 0, 0, 0, ...  ai+2
    si+5 = ...0, 0, 0, 0, h3,g3,h1,g1, 0, 0, 0, ...  ci+2
    si+6 = ...0, 0, 0, 0, 0, 0, h2,g2,h0,g0, 0, ...  ai+3
    si+7 = ...0, 0, 0, 0, 0, 0, h3,g3,h1,g1, 0, ...  ci+3
  

Using a standard dot product is grossly inefficient since most of the operands are zero. In practice the wavelet coefficient values are moved along the signal vector and a four element dot product is calculated. Expressed in terms of arrays, for the forward transform this would be:

  ai = s[i]*h0 + s[i+1]*h1 + s[i+2]*h2 + s[i+3]*h3
  ci = s[i]*g0 + s[i+1]*g1 + s[i+2]*g2 + s[i+3]*g3
  

This works fine if we have an infinite data set, since we don't have to worry about shifting the coefficients "off the end" of the signal.

I sometimes joke that I left my infinite data set in my other bear suit. The only problem with the algorithm described so far is that we don't have an infinite signal. The signal is finite. In fact not only must the signal be finite, but it must have a power of two number of elements.

If i=N-1, the i+2 and i+3 elements will be beyond the end of the array. There are a number of methods for handling the wavelet edge problem. This version of the algorithm acts like the data is periodic, where the data at the start of the signal wraps around to the end.

This algorithm uses a temporary array. A Lifting Scheme version of the Daubechies D4 algorithm does not require a temporary. The matrix discussion above is based on material from Ripples in Mathematics, by Jensen and Cour-Harbo. Any error are mine.

Author: Ian Kaplan
Use: You may use this software for any purpose as long as I cannot be held liable for the result. Please credit me with authorship if use use this source code.

This comment is formatted for the doxygen documentation generator

Definition at line 132 of file daub.h.


Constructor & Destructor Documentation

template<class T>
Daubechies<T>::Daubechies<T> ( ) [inline]
 

Initialize the filter constants used in the Daubechies D4 transform.

Definition at line 243 of file daub.h.

00244   {
00245     const double sqrt_3 = sqrt( 3 );
00246     const double denom = 4 * sqrt( 2 );
00247     
00248     //
00249     // forward transform scaling (smoothing) coefficients
00250     //
00251     h0 = (1 + sqrt_3)/denom;
00252     h1 = (3 + sqrt_3)/denom;
00253     h2 = (3 - sqrt_3)/denom;
00254     h3 = (1 - sqrt_3)/denom;
00255     //
00256     // forward transform wavelet coefficients (a.k.a. high
00257     // pass filter coefficients)
00258     //
00259     g0 =  h3;
00260     g1 = -h2;
00261     g2 =  h1;
00262     g3 = -h0;
00263     
00264     Ih0 = h2;
00265     Ih1 = g2;  // h1
00266     Ih2 = h0;
00267     Ih3 = g0;  // h3
00268     
00269     Ig0 = h3;
00270     Ig1 = g3;  // -h0
00271     Ig2 = h1;
00272     Ig3 = g1;  // -h2
00273   }


Member Function Documentation

template<class T>
void Daubechies<T>::forwardStep ( T & a,
const int n ) [inline, virtual]
 

Forward Daubechies D4 transform step.

Reimplemented from liftbase.

Definition at line 160 of file daub.h.

00161   {
00162     if (n >= 4) {
00163       int i, j;
00164       const int half = n >> 1;
00165       
00166       double* tmp = new double[n];
00167       
00168       for (i = 0, j = 0; j < n-3; j += 2, i++) {
00169         tmp[i]      = a[j]*h0 + a[j+1]*h1 + a[j+2]*h2 + a[j+3]*h3;
00170         tmp[i+half] = a[j]*g0 + a[j+1]*g1 + a[j+2]*g2 + a[j+3]*g3;
00171       }
00172       
00173       tmp[i]      = a[n-2]*h0 + a[n-1]*h1 + a[0]*h2 + a[1]*h3;
00174       tmp[i+half] = a[n-2]*g0 + a[n-1]*g1 + a[0]*g2 + a[1]*g3;
00175       
00176       for (i = 0; i < n; i++) {
00177         a[i] = tmp[i];
00178       }
00179       delete [] tmp;
00180     }
00181   }

template<class T>
void Daubechies<T>::forwardStepRev ( T & a,
const int n ) [inline, virtual]
 

Forward Daubechies D4 transform step, where the locations for the high and low pass filters are reversed.

Reimplemented from liftbase.

Definition at line 188 of file daub.h.

00189   {
00190     if (n >= 4) {
00191       int i, j;
00192       const int half = n >> 1;
00193       
00194       double* tmp = new double[n];
00195       
00196       for (i = 0, j = 0; j < n-3; j += 2, i++) {
00197         tmp[i+half] = a[j]*h0 + a[j+1]*h1 + a[j+2]*h2 + a[j+3]*h3;
00198         tmp[i]      = a[j]*g0 + a[j+1]*g1 + a[j+2]*g2 + a[j+3]*g3;
00199       }
00200       
00201       tmp[i+half] = a[n-2]*h0 + a[n-1]*h1 + a[0]*h2 + a[1]*h3;
00202       tmp[i]      = a[n-2]*g0 + a[n-1]*g1 + a[0]*g2 + a[1]*g3;
00203       
00204       for (i = 0; i < n; i++) {
00205         a[i] = tmp[i];
00206       }
00207       delete [] tmp;
00208     }
00209   }

template<class T>
void Daubechies<T>::forwardTrans ( T & ts,
int N ) [inline, virtual]
 

Forward Daubechies D4 transform.

Reimplemented from liftbase.

Definition at line 278 of file daub.h.

00279   {
00280     int n;
00281     for (n = N; n >= 4; n >>= 1) {
00282       forwardStep( ts, n );
00283     }
00284   } // forwardTrans

template<class T>
void Daubechies<T>::inverseStep ( T & a,
const int n ) [inline, virtual]
 

Inverse Daubechies D4 transform.

Reimplemented from liftbase.

Definition at line 214 of file daub.h.

00215   {
00216     if (n >= 4) {
00217       int i, j;
00218       const int half = n >> 1;
00219       const int halfPls1 = half + 1;
00220       
00221       double* tmp = new double[n];
00222       
00223       //      last smooth val  last coef.  first smooth  first coef
00224       tmp[0] = a[half-1]*Ih0 + a[n-1]*Ih1 + a[0]*Ih2 + a[half]*Ih3;
00225       tmp[1] = a[half-1]*Ig0 + a[n-1]*Ig1 + a[0]*Ig2 + a[half]*Ig3;
00226       for (i = 0, j = 2; i < half-1; i++) {
00227         //     smooth val     coef. val       smooth val    coef. val
00228         tmp[j++] = a[i]*Ih0 + a[i+half]*Ih1 + a[i+1]*Ih2 + a[i+halfPls1]*Ih3;
00229         tmp[j++] = a[i]*Ig0 + a[i+half]*Ig1 + a[i+1]*Ig2 + a[i+halfPls1]*Ig3;
00230       }
00231       for (i = 0; i < n; i++) {
00232         a[i] = tmp[i];
00233       }
00234       delete [] tmp;
00235     }
00236   } // inverseStep

template<class T>
void Daubechies<T>::inverseTrans ( T & coef,
int N ) [inline, virtual]
 

Inverse Daubechies D4 transform.

Reimplemented from liftbase.

Definition at line 290 of file daub.h.

00291   {
00292     int n;
00293     for (n = 4; n <= N; n <<= 1) {
00294       inverseStep( coef, n );
00295     }
00296   } // inverseTrans

template<class T>
void Daubechies<T>::predict ( T & vec,
int N,
transDirection direction ) [inline, protected, virtual]
 

Predict step, to be defined by the subclass.

Parameters:
vec   input array
N   size of region to act on (from 0..N-1)
direction   forward or inverse transform

Reimplemented from liftbase.

Definition at line 135 of file daub.h.

00136   {
00137     assert( false );
00138   } // predict

template<class T>
void Daubechies<T>::update ( T & vec,
int N,
transDirection direction ) [inline, protected, virtual]
 

Update step, to be defined by the subclass.

Parameters:
vec   input array
N   size of region to act on (from 0..N-1)
direction   forward or inverse transform

Reimplemented from liftbase.

Definition at line 140 of file daub.h.

00141   {
00142     assert( false );
00143   } // update


Member Data Documentation

template<class T>
double Daubechies<T>::Ig0 [private]
 

Definition at line 152 of file daub.h.

template<class T>
double Daubechies<T>::Ig1 [private]
 

Definition at line 152 of file daub.h.

template<class T>
double Daubechies<T>::Ig2 [private]
 

Definition at line 152 of file daub.h.

template<class T>
double Daubechies<T>::Ig3 [private]
 

Definition at line 152 of file daub.h.

template<class T>
double Daubechies<T>::Ih0 [private]
 

Definition at line 151 of file daub.h.

template<class T>
double Daubechies<T>::Ih1 [private]
 

Definition at line 151 of file daub.h.

template<class T>
double Daubechies<T>::Ih2 [private]
 

Definition at line 151 of file daub.h.

template<class T>
double Daubechies<T>::Ih3 [private]
 

Definition at line 151 of file daub.h.

template<class T>
double Daubechies<T>::g0 [private]
 

forward transform wave coefficients.

Definition at line 149 of file daub.h.

template<class T>
double Daubechies<T>::g1 [private]
 

forward transform wave coefficients.

Definition at line 149 of file daub.h.

template<class T>
double Daubechies<T>::g2 [private]
 

forward transform wave coefficients.

Definition at line 149 of file daub.h.

template<class T>
double Daubechies<T>::g3 [private]
 

forward transform wave coefficients.

Definition at line 149 of file daub.h.

template<class T>
double Daubechies<T>::h0 [private]
 

forward transform scaling coefficients.

Definition at line 147 of file daub.h.

template<class T>
double Daubechies<T>::h1 [private]
 

forward transform scaling coefficients.

Definition at line 147 of file daub.h.

template<class T>
double Daubechies<T>::h2 [private]
 

forward transform scaling coefficients.

Definition at line 147 of file daub.h.

template<class T>
double Daubechies<T>::h3 [private]
 

forward transform scaling coefficients.

Definition at line 147 of file daub.h.


The documentation for this class was generated from the following file:
Generated at Thu May 22 21:12:35 2003 for Hurst Exponent Calculation and Supporting Statistics by doxygen1.2.8.1 written by Dimitri van Heesch, © 1997-2001