Main Page   Compound List   File List   Compound Members  

Daubechies Class Reference

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

#include <daub.h>

List of all members.

Public Methods

 Daubechies ()
void daubTrans (double *ts, int N)
void invDaubTrans (double *coef, int N)

Private Methods

void transform (double *a, const int n)
void invTransform (double *a, const int n)

Private Attributes

double h0
double h1
double h2
double h3
double g0
double g1
double g2
double g3
double Ih0
double Ih1
double Ih2
double Ih3
double Ig0
double Ig1
double Ig2
double Ig3


Detailed Description

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


Constructor & Destructor Documentation

Daubechies::Daubechies ( ) [inline]
 

00215    {
00216       const double sqrt_3 = sqrt( 3 );
00217       const double denom = 4 * sqrt( 2 );
00218 
00219       //
00220       // forward transform scaling (smoothing) coefficients
00221       //
00222       h0 = (1 + sqrt_3)/denom;
00223       h1 = (3 + sqrt_3)/denom;
00224       h2 = (3 - sqrt_3)/denom;
00225       h3 = (1 - sqrt_3)/denom;
00226       //
00227       // forward transform wavelet coefficients
00228       //
00229       g0 =  h3;
00230       g1 = -h2;
00231       g2 =  h1;
00232       g3 = -h0;
00233 
00234       Ih0 = h2;
00235       Ih1 = g2;  // h1
00236       Ih2 = h0;
00237       Ih3 = g0;  // h3
00238 
00239       Ig0 = h3;
00240       Ig1 = g3;  // -h0
00241       Ig2 = h1;
00242       Ig3 = g1;  // -h2
00243    }


Member Function Documentation

void Daubechies::daubTrans ( double * ts,
int N ) [inline]
 

00246    {
00247       int n;
00248       for (n = N; n >= 4; n >>= 1) {
00249          transform( ts, n );
00250       }
00251    }

void Daubechies::invDaubTrans ( double * coef,
int N ) [inline]
 

00255    {
00256       int n;
00257       for (n = 4; n <= N; n <<= 1) {
00258          invTransform( coef, n );
00259       }
00260    }

void Daubechies::invTransform ( double * a,
const int n ) [inline, private]
 

Inverse Daubechies D4 transform.

00188    {
00189      if (n >= 4) {
00190        int i, j;
00191        const int half = n >> 1;
00192        const int halfPls1 = half + 1;
00193 
00194        double* tmp = new double[n];
00195 
00196        //      last smooth val  last coef.  first smooth  first coef
00197        tmp[0] = a[half-1]*Ih0 + a[n-1]*Ih1 + a[0]*Ih2 + a[half]*Ih3;
00198        tmp[1] = a[half-1]*Ig0 + a[n-1]*Ig1 + a[0]*Ig2 + a[half]*Ig3;
00199        for (i = 0, j = 2; i < half-1; i++) {
00200          //     smooth val     coef. val       smooth val    coef. val
00201          tmp[j++] = a[i]*Ih0 + a[i+half]*Ih1 + a[i+1]*Ih2 + a[i+halfPls1]*Ih3;
00202          tmp[j++] = a[i]*Ig0 + a[i+half]*Ig1 + a[i+1]*Ig2 + a[i+halfPls1]*Ig3;
00203        }
00204        for (i = 0; i < n; i++) {
00205          a[i] = tmp[i];
00206        }
00207        delete [] tmp;
00208      }
00209    }

void Daubechies::transform ( double * a,
const int n ) [inline, private]
 

Forward Daubechies D4 transform.

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


Member Data Documentation

double Daubechies::Ig0 [private]
 

double Daubechies::Ig1 [private]
 

double Daubechies::Ig2 [private]
 

double Daubechies::Ig3 [private]
 

double Daubechies::Ih0 [private]
 

double Daubechies::Ih1 [private]
 

double Daubechies::Ih2 [private]
 

double Daubechies::Ih3 [private]
 

double Daubechies::g0 [private]
 

forward transform wave coefficients.

double Daubechies::g1 [private]
 

forward transform wave coefficients.

double Daubechies::g2 [private]
 

forward transform wave coefficients.

double Daubechies::g3 [private]
 

forward transform wave coefficients.

double Daubechies::h0 [private]
 

forward transform scaling coefficients.

double Daubechies::h1 [private]
 

forward transform scaling coefficients.

double Daubechies::h2 [private]
 

forward transform scaling coefficients.

double Daubechies::h3 [private]
 

forward transform scaling coefficients.


The documentation for this class was generated from the following file:
Generated at Wed Jan 23 20:56:40 2002 for Daubechies D4 Wavelet Transform by doxygen1.2.8.1 written by Dimitri van Heesch, © 1997-2001