/**
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
*/
template
class Daubechies : public liftbase {
protected:
void predict( T& vec, int N, transDirection direction )
{
assert( false );
} // predict
virtual void update( T& vec, int N, transDirection direction )
{
assert( false );
} // update
private:
/** forward transform scaling coefficients */
double h0, h1, h2, h3;
/** forward transform wave coefficients */
double g0, g1, g2, g3;
double Ih0, Ih1, Ih2, Ih3;
double Ig0, Ig1, Ig2, Ig3;
public:
/**
Forward Daubechies D4 transform step
*/
void forwardStep( T& a, const int n )
{
if (n >= 4) {
int i, j;
const int half = n >> 1;
double* tmp = new double[n];
for (i = 0, j = 0; j < n-3; j += 2, i++) {
tmp[i] = a[j]*h0 + a[j+1]*h1 + a[j+2]*h2 + a[j+3]*h3;
tmp[i+half] = a[j]*g0 + a[j+1]*g1 + a[j+2]*g2 + a[j+3]*g3;
}
tmp[i] = a[n-2]*h0 + a[n-1]*h1 + a[0]*h2 + a[1]*h3;
tmp[i+half] = a[n-2]*g0 + a[n-1]*g1 + a[0]*g2 + a[1]*g3;
for (i = 0; i < n; i++) {
a[i] = tmp[i];
}
delete [] tmp;
}
}
/**
Forward Daubechies D4 transform step, where the locations
for the high and low pass filters are reversed.
*/
void forwardStepRev( T& a, const int n )
{
if (n >= 4) {
int i, j;
const int half = n >> 1;
double* tmp = new double[n];
for (i = 0, j = 0; j < n-3; j += 2, i++) {
tmp[i+half] = a[j]*h0 + a[j+1]*h1 + a[j+2]*h2 + a[j+3]*h3;
tmp[i] = a[j]*g0 + a[j+1]*g1 + a[j+2]*g2 + a[j+3]*g3;
}
tmp[i+half] = a[n-2]*h0 + a[n-1]*h1 + a[0]*h2 + a[1]*h3;
tmp[i] = a[n-2]*g0 + a[n-1]*g1 + a[0]*g2 + a[1]*g3;
for (i = 0; i < n; i++) {
a[i] = tmp[i];
}
delete [] tmp;
}
}
/**
Inverse Daubechies D4 transform
*/
void inverseStep( T& a, const int n )
{
if (n >= 4) {
int i, j;
const int half = n >> 1;
const int halfPls1 = half + 1;
double* tmp = new double[n];
// last smooth val last coef. first smooth first coef
tmp[0] = a[half-1]*Ih0 + a[n-1]*Ih1 + a[0]*Ih2 + a[half]*Ih3;
tmp[1] = a[half-1]*Ig0 + a[n-1]*Ig1 + a[0]*Ig2 + a[half]*Ig3;
for (i = 0, j = 2; i < half-1; i++) {
// smooth val coef. val smooth val coef. val
tmp[j++] = a[i]*Ih0 + a[i+half]*Ih1 + a[i+1]*Ih2 + a[i+halfPls1]*Ih3;
tmp[j++] = a[i]*Ig0 + a[i+half]*Ig1 + a[i+1]*Ig2 + a[i+halfPls1]*Ig3;
}
for (i = 0; i < n; i++) {
a[i] = tmp[i];
}
delete [] tmp;
}
} // inverseStep
/**
Initialize the filter constants used in the Daubechies
D4 transform.
*/
Daubechies()
{
const double sqrt_3 = sqrt( 3 );
const double denom = 4 * sqrt( 2 );
//
// forward transform scaling (smoothing) coefficients
//
h0 = (1 + sqrt_3)/denom;
h1 = (3 + sqrt_3)/denom;
h2 = (3 - sqrt_3)/denom;
h3 = (1 - sqrt_3)/denom;
//
// forward transform wavelet coefficients (a.k.a. high
// pass filter coefficients)
//
g0 = h3;
g1 = -h2;
g2 = h1;
g3 = -h0;
Ih0 = h2;
Ih1 = g2; // h1
Ih2 = h0;
Ih3 = g0; // h3
Ig0 = h3;
Ig1 = g3; // -h0
Ig2 = h1;
Ig3 = g1; // -h2
}
/**
Forward Daubechies D4 transform
*/
void forwardTrans( T& ts, int N )
{
int n;
for (n = N; n >= 4; n >>= 1) {
forwardStep( ts, n );
}
} // forwardTrans
/**
Inverse Daubechies D4 transform
*/
void inverseTrans( T& coef, int N )
{
int n;
for (n = 4; n <= N; n <<= 1) {
inverseStep( coef, n );
}
} // inverseTrans
}; // Daubechies
packet/include/fifo_list.h 0100644 0000764 0000764 00000011177 07467262147 014564 0 ustar iank iank
#ifndef FIFO_LIST_H
#define FIFO_LIST_H
/** \file
The documentation in this file is formatted for doxygen
(see www.doxygen.org).
Copyright and Use
You may use this source code without limitation and without
fee as long as you include:
This software was written and is copyrighted by Ian Kaplan, Bear
Products International, www.bearcave.com, 2002.
This software is provided "as is", without any warranty or
claim as to its usefulness. Anyone who uses this source code
uses it at their own risk. Nor is any support provided by
Ian Kaplan and Bear Products International.
Please send any bug fixes or suggested source changes to:
iank@bearcave.com
@author Ian Kaplan
*/
#include "blockpool.h"
/**
template class FIFO_LIST
This is a generic list type for a list that has both a head and a
tail pointer. In this list, items are added to the tail. When
read from the front of the list, items will be read in a first-in,
first-out order (FIFO). The template should be instantiated with a
scalar type, like an integer or a pointer to a larger type (e.g., a
string or a structure). For example
FIFO_LIST list;
FIFO_LIST list;
FIFO_LIST list;
The list "links" are allocated in a memory pool. This allocation
is handled by the local version of new in list_type. The memory
allocated for the list_type objects will be deallocated when
the memory pool is freed.
*/
template
class FIFO_LIST
{
public:
/**
List backbone class
*/
class list_type {
public:
/** data element */
T data;
/** pointer to the next element */
list_type *next;
/** override the default new operator to allocate
list_type objects from the a memory pool */
void *operator new(size_t num_bytes)
{
block_pool mem_pool;
void *mem_addr = mem_pool.pool_alloc( num_bytes );
return mem_addr;
} // new
}; // class list_type
private:
/** list head */
list_type *list;
/** list tail (items are added to the tail) */
list_type *tail;
public:
/** define a handle type to abstract the list_type type */
typedef list_type *handle;
public:
/** class constructor */
FIFO_LIST(void)
{
list = 0;
tail = 0;
}
/** default destructor does nothing */
~FIFO_LIST(void) {}
/** deallocate the list */
void dealloc(void)
{
while ( remove() != 0 )
/* nada */;
} // dealloc
/** add an element to the FIFO list */
void add( T data )
{
list_type *t;
t = new list_type();
t->data = data;
t->next = 0;
if (list == 0) {
list = t;
tail = t;
}
else {
tail->next = t;
tail = t;
}
} // add
/** reverse the list */
void reverse(void)
{
list_type *elem, *prev, *next;
prev = 0;
next = 0;
tail = list;
for (elem = list; elem != 0; prev = elem, elem = next) {
next = elem->next;
elem->next = prev;
} // for
list = prev;
} // reverse
/** return the lenght of the list */
unsigned int length(void)
{
list_type *elem;
unsigned int cnt = 0;
for (elem = list; elem != 0; elem = elem->next)
cnt++;
return cnt;
} // lenght
/**
remove
Remove an element from the start of the list and return the first
element of the remaining list.
This function relies on the fact that the list elements were
allocated from a pool. The memory for these elements will be
recovered when the pool is deallocated.
*/
handle remove(void)
{
list_type *t;
if (list != 0) {
t = list;
list = t->next;
// no delete t;
}
if (list == 0)
tail = 0;
return list;
} // remove
/** given a handle, return the associated data item */
T get_item( handle h)
{
return h->data;
} // get_item
/** get the first element from the list */
handle first(void)
{
return list;
} // first
/** return the last element in the list */
handle last(void)
{
return tail;
} // last
/** iterator to get the next element */
handle next(handle h)
{
list_type *next = 0;
if (h != 0) {
next = h->next;
}
return next;
} // next
}; // template class FIFO_LIST
#endif
packet/include/grow_array.h 0100644 0000764 0000764 00000011125 07525053166 014745 0 ustar iank iank
#ifndef _GROW_ARRAY_H_
#define _GROW_ARRAY_H_
/** \file
The documentation in this file is formatted for doxygen
(see www.doxyeng.org).
Copyright and Use
You may use this source code without limitation and without
fee as long as you include:
This software was written and is copyrighted by Ian Kaplan, Bear
Products International, www.bearcave.com, 2002.
This software is provided "as is", without any warranty or
claim as to its usefulness. Anyone who uses this source code
uses it at their own risk. Nor is any support provided by
Ian Kaplan and Bear Products International.
Please send any bug fixes or suggested source changes to:
iank@bearcave.com
@author Ian Kaplan
*/
#include "blockpool.h"
/**
This file defines an array template class that will grow as
elements are added to the end of the array. This is similar
to the STL class.
This array class is designed for dense arrays where the all
elements of the array are used.
Usage:
- Elements are added to the end of the array via the append
function.
- Elements in the array can be accessed via the [] operator.
- If the elements in the array are dynamicly allocated, the
user is responsible for deallocating these elements.
A doubling algorithm is used when the data size is expanded because
it minimizes the amount of copying that must be done. The array
will quickly grow to a the size needed to accomodate the data set
and no more copying will be necessary. For large arrays there is the
drawback that more memory may be allocated than is needed, since
the amount of memory used grows exponentially.
\author Ian Kaplan
*/
template
class GrowableArray {
private:
typedef enum { StartArraySize = 128 } bogus;
/** number of data elements */
size_t num_elem;
/** Array size (always <= num_elem) */
size_t array_size;
T *pArray;
private:
/**
twice
Double the amount of memory allocated for the array.
Return true if memory allocation succeeded, false
otherwise.
*/
bool twice()
{
bool rslt;
T *old_array = pArray;
size_t new_size = array_size * 2;
pArray = new T [ new_size ];
if (pArray != 0) {
rslt = true;
for (int i = 0; i < array_size; i++) {
pArray[i] = old_array[i];
}
delete [] old_array;
array_size = new_size;
}
else {
rslt = false;
}
return rslt;
} // twice
public:
GrowableArray()
{
pArray = new T[ StartArraySize ];
num_elem = 0;
array_size = StartArraySize;
} // GrowableArray constructor
/** destructor */
~GrowableArray()
{
if (pArray != NULL) {
delete [] pArray;
}
} // GrowableArray destructor
/** Length of the data (which is not necessarily the same
as the length of the internal array */
const size_t length(void) const { return num_elem; }
/** set array to zero length */
void set_to_zero()
{
num_elem = 0;
}
/** LHS [] operator */
T &operator[](const size_t i)
{
assert( i < num_elem );
return pArray[ i ];
}
/** RHS [] operator */
T operator[](const size_t i ) const
{
assert( i < num_elem );
return pArray[ i ];
}
/** Get a pointer to the internal data array */
const T *getData() const { return pArray; }
/** append an item to the end of the array */
void append( T item )
{
if (num_elem == array_size) {
bool allocOK = twice();
assert( allocOK );
}
pArray[ num_elem ] = item;
num_elem++;
} // append
/**
expand
Expand the number of array data slots by "amount" elements.
Note that "array_size" is the total amount of storage available
for data slots. "num_elem" is the number of data slots.
The bounds over which the array can be indexed is governed
by num_elem. Note that after expand() is called the new
data elements can be read, but their value is undefined until
they are initialized.
*/
void expand( size_t amount )
{
bool allocOK = true;
while (allocOK && num_elem + amount >= array_size) {
allocOK = twice();
assert( allocOK );
}
num_elem += amount;
} // expand
/** Remove one item from the end of the array. */
void remove(void)
{
if (num_elem > 0)
num_elem--;
}
}; // GrowableArray
#endif
packet/include/haar.h 0100644 0000764 0000764 00000014557 07525053205 013512 0 ustar iank iank
#ifndef _HAAR_H_
#define _HAAR_H_
#include
#include "liftbase.h"
/** \file
The documentation in this file is formatted for doxygen
(see www.doxygen.org).
Copyright and Use
You may use this source code without limitation and without
fee as long as you include:
This software was written and is copyrighted by Ian Kaplan, Bear
Products International, www.bearcave.com, 2002.
This software is provided "as is", without any warranty or
claim as to its usefulness. Anyone who uses this source code
uses it at their own risk. Nor is any support provided by
Ian Kaplan and Bear Products International.
Please send any bug fixes or suggested source changes to:
iank@bearcave.com
@author Ian Kaplan
*/
/**
Haar (flat line) wavelet.
As with all Lifting scheme wavelet transform functions, the
first stage of a transform step is the split stage. The
split step moves the even element to the first half of an
N element region and the odd elements to the second half of the N
element region.
The Lifting Scheme version of the Haar transform uses a wavelet
function (predict stage) that "predicts" that an odd element will
have the same value as it preceeding even element. Stated another
way, the odd element is "predicted" to be on a flat (zero slope
line) shared with the even point. The difference between this
"prediction" and the actual odd value replaces the odd element.
The wavelet scaling function (a.k.a. smoothing function) used
in the update stage calculates the average between an even and
an odd element.
The merge stage at the end of the inverse transform interleaves
odd and even elements from the two halves of the array
(e.g., ordering them even0, odd0,
even1, odd1, ...)
This is a template version of the Haar wavelet. The template must
be instantiated with an array or an object that acts like an array.
Objects that act like arrays define the left hand side and right
hand side index operators: [].
See www.bearcave.com for more information on wavelets and the
wavelet lifting scheme.
\author Ian Kaplan
*/
template
class haar : public liftbase {
protected:
/**
Haar predict step
*/
void predict( T& vec, int N, transDirection direction )
{
int half = N >> 1;
for (int i = 0; i < half; i++) {
double predictVal = vec[i];
int j = i + half;
if (direction == forward) {
vec[j] = vec[j] - predictVal;
}
else if (direction == inverse) {
vec[j] = vec[j] + predictVal;
}
else {
printf("haar::predict: bad direction value\n");
}
}
}
/**
Update step of the Haar wavelet transform.
The wavelet transform calculates a set of detail or
difference coefficients in the predict step. These
are stored in the upper half of the array. The update step
calculates an average from the even-odd element pairs.
The averages will replace the even elements in the
lower half of the array.
The Haar wavelet calculation used in the Lifting Scheme
is
dj+1, i = oddj+1, i = oddj, i - evenj, i
aj+1, i = evenj, i = (evenj, i + oddj, i)/2
Note that the Lifting Scheme uses an in-place algorithm. The odd
elements have been replaced by the detail coefficients in the
predict step. With a little algebra we can substitute the
coefficient calculation into the average calculation, which
gives us
aj+1, i = evenj, i = evenj, i + (oddj, i/2)
*/
void update( T& vec, int N, transDirection direction )
{
int half = N >> 1;
for (int i = 0; i < half; i++) {
int j = i + half;
double updateVal = vec[j] / 2.0;
if (direction == forward) {
vec[i] = vec[i] + updateVal;
}
else if (direction == inverse) {
vec[i] = vec[i] - updateVal;
}
else {
printf("update: bad direction value\n");
}
}
}
/**
The normalization step assures that each step of the wavelet
transform has the constant "energy" where energy is defined as
double energy = 0.0;
for (int n = 0; n < N; n++) {
energy = energy + (a[i] * a[i]);
}
See 5.2.1 of Ripples in Mathematics by Jensen
and la Cour-Harbo
The most common implementation of the Haar transform leaves out
the normalization step, since it does not make much of a
difference in many cases. However, in the case of the wavelet
packet transform, many of the cost functions are squares, so
normalization produces smaller wavelet values (although the
scaling function values are larger). This may lead to a better
wavelet packet result (e.g., a few large values and lots of small
values).
Normalization does have the disadvantage of destroying the
averaging property of the Haar wavelet algorithm. That is, the
final scaling factor is no longer the mean of the time series.
*/
void normalize( T& vec, int N, transDirection direction )
{
const double sqrt2 = sqrt( 2.0 );
int half = N >> 1;
for (int i = 0; i < half; i++) {
int j = i + half;
if (direction == forward) {
vec[i] = sqrt2 * vec[i];
vec[j] = vec[j]/sqrt2;
}
else if (direction == inverse) {
vec[i] = vec[i]/sqrt2;
vec[j] = sqrt2 * vec[j];
}
else {
printf("normalize: bad direction value\n");
}
} // for
} // normalize
/**
One inverse wavelet transform step, with normalization
*/
void inverseStep( T& vec, const int n )
{
normalize( vec, n, inverse );
update( vec, n, inverse );
predict( vec, n, inverse );
merge( vec, n );
} // inverseStep
/**
One step in the forward wavelet transform, with normalization
*/
void forwardStep( T& vec, const int n )
{
split( vec, n );
predict( vec, n, forward );
update( vec, n, forward );
normalize( vec, n, forward );
} // forwardStep
}; // haar
#endif
packet/include/haar_classic.h 0100644 0000764 0000764 00000010545 07525053227 015210 0 ustar iank iank
#ifndef _HAAR_CLASSIC_H_
#define _HAAR_CLASSIC_H_
#include
#include "liftbase.h"
/** \file
The documentation in this file is formatted for doxygen
(see www.doxygen.org)
Copyright and Use
You may use this source code without limitation and without
fee as long as you include:
This software was written and is copyrighted by Ian Kaplan, Bear
Products International, www.bearcave.com, 2001.
This software is provided "as is", without any warranty or
claim as to its usefulness. Anyone who uses this source code
uses it at their own risk. Nor is any support provided by
Ian Kaplan and Bear Products International.
Please send any bug fixes or suggested source changes to:
iank@bearcave.com
@author Ian Kaplan
*/
/**
A version of the classic Haar wavelet transform
This particular version of the Haar wavelet transform is frequently
given as the definition for the Haar transform. This version differs
from the lifting scheme version. In the case of the lifting
scheme version of the Haar transform, the inverse transform is
a mirror of the forward transform. The only difference is that
the plus and minus operators are interchanged. This algorithm
does not have this symmetry.
For a data set of N elements, a wavelet transform will calculate
N/2 smoothed values and N/2 difference values. In wavelet
terminology the smoothed values are calculated by the scaling
function and the difference (or coefficient) values are calculated
by the wavelet function.
This class implements one version of the Haar wavelet transform.
This particular version is used in the chapter 8 of Ripples
in Mathematics by Jensen and la Cour-Harbo to illustrate
the wavelet packet transform. I have used it to verify my version
of the wavelet packet algorithm.
In the description below, an element ai is an even element
and an element bi is an odd element.
In this version of the Haar wavelet transform the scaling (or
smoothing) function is
s = (a + b)/2
The wavelet function is
d = (a - b)/2
A lifting scheme expression is used in this implementation. Here
the wavelet function is calculated first. The wavelet results
overwrite the odd bi values. This means that the
smoothing function values must be calculated with the result
of the wavelet function. To recover the value of bi
we use the expression
b = a - 2d
s = (a + (a - 2d))/2
s = (2a - 2d)/2
s = a - d
The lifting scheme terminology is maintained in the algorithm,
although it does not fully apply.
This is a template version of the Haar wavelet. The template must
be instantiated with an array or an object that acts like an array.
Objects that act like arrays define the left hand side and right
hand side index operators: [].
See www.bearcave.com for more information on wavelets and the
wavelet lifting scheme.
\author Ian Kaplan
*/
template
class haar_classic : public liftbase {
protected:
/**
Calculate the Haar wavelet or difference function (high
pass filter)
*/
void predict( T& vec, int N, transDirection direction )
{
int half = N >> 1;
int cnt = 0;
for (int i = 0; i < half; i++) {
double predictVal = vec[i];
int j = i + half;
if (direction == forward) {
vec[j] = (predictVal - vec[j] )/2;
}
else if (direction == inverse) {
vec[j] = predictVal - (2 * vec[j]);
}
else {
printf("haar_classic::predict: bad direction value\n");
}
}
} // predict
/**
Calculate the smoothing or scaling function (low pass
filter)
*/
void update( T& vec, int N, transDirection direction )
{
int half = N >> 1;
for (int i = 0; i < half; i++) {
int j = i + half;
double updateVal = vec[j];
if (direction == forward) {
vec[i] = vec[i] - updateVal;
}
else if (direction == inverse) {
vec[i] = vec[i] + updateVal;
}
else {
printf("update: bad direction value\n");
}
}
} // update
}; // haar_classic
#endif
packet/include/haar_classicFreq.h 0100644 0000764 0000764 00000011031 07525053246 016016 0 ustar iank iank
#ifndef _HAAR_CLASSICFREQ_H_
#define _HAAR_CLASSICFREQ_H_
#include "haar_classic.h"
/**
haar_classicFreq
An extension of the "Haar classic" algorithm for frequency
analysis.
The haar_classicFreq template extends the haar_classic
templace. The haar_classic template implements what
I call the "Haar classic" algorithm. Here the high
pass filter (wavelet) is
b'i = (ai - bi)/2
where ai is an even element and bi is it odd
element neighbor (of course using a lifting scheme style
implementation, the even elements are moved to the first half of the
array and the odd elements are moved to the second half.
In the Haar classic algorithm the high pass result is placed
in the upper half of the array. In the haar_classicFreq version
the result is placed in the lower half of the array and the
equation becomes
a'i = (ai - bi)/2
The "Haar classic" low pass filter (smoothing function) is
a'i = (ai + bi)/2
In the Haar classic algorithm the low pass filter result is placed
in the lower half of the array. In the haar_classicFreq version the
result is placed in the upper half of the array and the equation
becomes
b'i = (ai + bi)/2
The calculation of the high pass filter is done first and
overwrites the even elements (e.g., a'i). To
recover the ai values given a'i and
bi:
a'i = (ai - bi)/2
2 * a'i = ai - bi
ai = 2 * a'i + bi
Substituting this into the equation for calculating
b'i we get:
b'i = ((2 * a'i + bi) + bi)/2
b'i = (2 * a'i + 2 * bi)/2
b'i = a'i + bi
These equations differ from the Haar classic forward transform
equations.
\author Ian Kaplan
*/
template
class haar_classicFreq : public haar_classic {
protected:
/**
In the standard wavelet transform, the high pass filter is applied
to the upper half of the array (this is the predict phase, in
lifting scheme terminology).
In the reverse transform step, which is used for wavelet frequency
analysis, the high pass filter is applied to the lower half of the
array.
*/
void predictRev( T& vec, int N, transDirection direction )
{
int half = N >> 1;
int cnt = 0;
for (int i = 0; i < half; i++) {
int j = i + half;
if (direction == forward) {
vec[i] = (vec[i] - vec[j] )/2;
}
else if (direction == inverse) {
vec[i] = (2 * vec[i]) + vec[j];
}
else {
printf("predictRev: bad direction value\n");
}
}
} // predictRev
/**
Reverse low pass filter
In the standard wavelet transform the low pass filter is
applied to the data set (consisting of N elements) and
the result is placed in the lower half of the array
(the lower N/2 elements).
In the reverse transform step, which is used for wavelet
frequency analysis, the result of the low pass filter
is placed in the upper half of the array (the upper
N/2 elements).
*/
void updateRev( T& vec, int N, transDirection direction )
{
int half = N >> 1;
for (int i = 0; i < half; i++) {
int j = i + half;
if (direction == forward) {
vec[j] = vec[j] + vec[i];
}
else if (direction == inverse) {
vec[j] = vec[j] - vec[i];
}
else {
printf("updateRev: bad direction value\n");
}
}
} // updateRev
public:
/**
One forward step of the reverse Haar classic transform, where the
results for the high and low pass filters are reversed.
*/
void forwardStepRev( T& vec, const int n )
{
split( vec, n );
predictRev( vec, n, forward );
updateRev( vec, n, forward );
}
/**
One inverse step of the reverse Haar classic transform, where the
results for the high and low pass filters are reversed.
*/
void inverseStepRev( T& vec, const int n )
{
updateRev( vec, n, inverse );
predictRev( vec, n, inverse );
merge( vec, n );
}
}; // haar_classicFreq
#endif
packet/include/invpacktree.h 0100644 0000764 0000764 00000006326 07525053262 015110 0 ustar iank iank
#ifndef _INVPACKTREE_H_
#define _INVPACKTREE_H_
/** \file
The documentation in this file is formatted for doxygen
(see www.doxygen.org).
Copyright and Use
You may use this source code without limitation and without
fee as long as you include:
This software was written and is copyrighted by Ian Kaplan, Bear
Products International, www.bearcave.com, 2002.
This software is provided "as is", without any warranty or
claim as to its usefulness. Anyone who uses this source code
uses it at their own risk. Nor is any support provided by
Ian Kaplan and Bear Products International.
Please send any bug fixes or suggested source changes to:
iank@bearcave.com
@author Ian Kaplan
*/
#include "liftbase.h"
#include "list.h"
#include "packcontainer.h"
#include "packdata.h"
#include "packdata_list.h"
/**
Inverse wavelet packet transform
The invpacktree constructor is passed a packdata_list object and a
wavelet transform object. It calculates the inverse wavelet
packet transform, using the data in the packdata_list object
and the inverse wavelet transform step function of the wavelet
transform object. The best basis data is destroyed in calculating
the inverse transform.
The packdata_list object contains the "best basis" result from a
wavelet packet transform. The wavelet packet transform should have
been calculated with the same wavelet transform as the object passed
to this constructor.
After the constructor completes, the data result can be
obtained by calling the getData() function.
The wavelet transforms used by this object are all derived
from the liftbase class and are "lifting scheme" wavelet
transforms.
I have found the wavelet literature difficult, in general. When it
comes to the wavelet packet transform I have found most of it
impossible to understand. Impossible, that it until I got the book
Ripples in Mathematics by Jensen and la Cour-Harbo, Springer
Verlag, 2001. The wavelet packet transform, for which this class
is the inverse, is heavily based on Chapter 8 of Ripples in
Mathematics.
I have found very little material on the actual implementation of
the inverse wavelet packet transform. This algorithm is my own
design.
\author Ian Kaplan
*/
class invpacktree {
private:
/** wavelet transform object */
liftbase *waveObj;
/** disallow the copy constructor */
invpacktree( const invpacktree &rhs ) {}
/** inverse wavelet packet transform calculation stack */
LIST stack;
/** pointer to the inverse transform result */
const double *data;
/** length of data */
size_t N;
private:
void new_level( packdata *elem );
void add_elem( packdata *elem );
void reduce();
public:
invpacktree( packdata_list &list,
liftbase *w );
/** The destructor does nothing */
~invpacktree() {}
/** Get the result of the inverse packet transform */
const double *getData() { return data; }
void pr();
};
#endif
packet/include/liftbase.h 0100644 0000764 0000764 00000023074 07525053311 014360 0 ustar iank iank
#ifndef _LIFTBASE_H_
#define _LIFTBASE_H_
/** \file
Copyright and Use
You may use this source code without limitation and without
fee as long as you include:
This software was written and is copyrighted by Ian Kaplan, Bear
Products International, www.bearcave.com, 2002.
This software is provided "as is", without any warranty or
claim as to its usefulness. Anyone who uses this source code
uses it at their own risk. Nor is any support provided by
Ian Kaplan and Bear Products International.
Please send any bug fixes or suggested source changes to:
iank@bearcave.com
@author Ian Kaplan
*/
#include
/**
This is the base class for simple Lifting Scheme wavelets using
split, predict, update or update, predict, merge steps.
Simple lifting scheme wavelets consist of three steps,
a split/merge step, predict step and an update step:
-
The split step divides the elements in an array so that
the even elements are in the first half and the odd
elements are in the second half.
-
The merge step is the inverse of the split step. It takes
two regions of an array, an odd region and an even region
and merges them into a new region where an even element
alternates with an odd element.
-
The predict step calculates the difference
between an odd element and its predicted value based
on the even elements. The difference between the
predicted value and the actual value replaces the
odd element.
-
The predict step operates on the odd elements. The update
step operates on the even element, replacing them with a
difference between the predict value and the actual odd element.
The update step replaces each even element with an average.
The result of the update step becomes the input to the
next recursive step in the wavelet calculation.
The split and merge methods are shared by all Lifting Scheme wavelet
algorithms. This base class provides the transform and inverse
transform methods (forwardTrans and inverseTrans). The predict
and update methods are abstract and are defined for a particular
Lifting Scheme wavelet sub-class.
This is a template version of the lifting scheme base class. The
template must be instantiated with an array or an object that acts
like an array. Objects that act like arrays define the left hand
side and right hand side index operators: []. To allow wavelet
transforms based on this base class to be used with the wavelet
packet transform, this class makes public both the forward
and inverse transforms (forwardTrans and inverseTrans) and
the forward and inverse transform steps (forwardStep and
inverseStep). These "step" functions are used to calculate
the wavelet packet transform.
Instantiating the Template
The liftbase template takes two type arguments:
-
The type of the array or '[]' operator indexable object.
-
The type of the data element.
The simplest example is a wavelet class derived from an instance of
the liftbase tempate which takes a double array and has a double
element type. This declaration is shown below:
class Haar : public liftbase
An object type can be used for the first template argument,
as long as the object supports the '[]' operator, which returns
an element whose type is defined by the second argument. In the
example below, the packcontainer '[]' operator returns a
double.
class Poly : public liftbase
References:
\author Ian Kaplan
*/
template
class liftbase {
protected:
typedef enum {
/** "enumeration" for forward wavelet transform */
forward = 1,
/** "enumeration" for inverse wavelet transform */
inverse = 2
} transDirection;
/**
Split the vec into even and odd elements,
where the even elements are in the first half
of the vector and the odd elements are in the
second half.
*/
void split( T& vec, int N )
{
int start = 1;
int end = N - 1;
while (start < end) {
for (int i = start; i < end; i = i + 2) {
T_elem tmp = vec[i];
vec[i] = vec[i+1];
vec[i+1] = tmp;
}
start = start + 1;
end = end - 1;
}
}
/**
Merge the odd elements from the second half of the N element
region in the array with the even elements in the first
half of the N element region. The result will be the
combination of the odd and even elements in a region
of length N.
*/
void merge( T& vec, int N )
{
int half = N >> 1;
int start = half-1;
int end = half;
while (start > 0) {
for (int i = start; i < end; i = i + 2) {
T_elem tmp = vec[i];
vec[i] = vec[i+1];
vec[i+1] = tmp;
}
start = start - 1;
end = end + 1;
}
}
/**
Predict step, to be defined by the subclass
@param vec input array
@param N size of region to act on (from 0..N-1)
@param direction forward or inverse transform
*/
virtual void predict( T& vec, int N, transDirection direction ) = 0;
/**
Reverse predict step.
The predict step applied the high pass filter to the data
set and places the result in the upper half of the array.
The reverse predict step applies the high pass filter and
places the result in the lower half of the array.
This reverse predict step is only used by wavelet packet
frequency analysis algorithms. The default version
of this algorihtm does nothing.
*/
virtual void predictRev( T& vec, int N, transDirection direction ) {};
/**
Update step, to be defined by the subclass
@param vec input array
@param N size of region to act on (from 0..N-1)
@param direction forward or inverse transform
*/
virtual void update( T& vec, int N, transDirection direction ) = 0;
/**
Reverse update step
*/
virtual void updateRev( T& vec, int N, transDirection direction ) {}
public:
/**
One step in the forward wavelet transform
*/
virtual void forwardStep( T& vec, const int n )
{
split( vec, n );
predict( vec, n, forward );
update( vec, n, forward );
} // forwardStep
/**
Reverse forward transform step. The result of the high
pass filter is stored in the lower half of the array
and the result of the low pass filter is stored in the
upper half.
This function should be defined by any subclass that
is used for wavelet frequency analysis.
*/
virtual void forwardStepRev( T& vec, const int N )
{
assert(false);
}
/**
Simple wavelet Lifting Scheme forward transform
forwardTrans is passed an indexable object. The object must
contain a power of two number of data elements. Lifting Scheme
wavelet transforms are calculated in-place and the result is
returned in the argument array.
The result of forwardTrans is a set of wavelet coefficients
ordered by increasing frequency and an approximate average
of the input data set in vec[0]. The coefficient bands
follow this element in powers of two (e.g., 1, 2, 4, 8...).
*/
virtual void forwardTrans( T& vec, const int N )
{
for (int n = N; n > 1; n = n >> 1) {
forwardStep( vec, n );
}
} // forwardTrans
/**
One inverse wavelet transform step
*/
virtual void inverseStep( T& vec, const int n )
{
update( vec, n, inverse );
predict( vec, n, inverse );
merge( vec, n );
}
/**
Reverse inverse transform step. Calculate the inverse transform
from a high pass filter result stored in the lower half of the
array and a low pass filter result stored in the upper half.
This function should be defined by any subclass that
is used for wavelet frequency analysis.
*/
virtual void inverseStepRev( T& vec, const int n )
{
assert( false );
}
/**
Default two step Lifting Scheme inverse wavelet transform
inverseTrans is passed the result of an ordered wavelet
transform, consisting of an average and a set of wavelet
coefficients. The inverse transform is calculated
in-place and the result is returned in the argument array.
*/
virtual void inverseTrans( T& vec, const int N )
{
for (int n = 2; n <= N; n = n << 1) {
inverseStep( vec, n );
}
} // inverseTrans
}; // liftbase
#endif
packet/include/line.h 0100644 0000764 0000764 00000020037 07525053337 013522 0 ustar iank iank
#ifndef _LINE_H_
#define _LINE_H_
#include "liftbase.h"
/**
Line (with slope) wavelet
The wavelet Lifting Scheme "line" wavelet approximates the
data set using a line with with slope (in contrast to the
Haar wavelet where a line has zero slope is used to approximate
the data).
The predict stage of the line wavelet "predicts" that an odd point
will lie midway between its two neighboring even points. That is,
that the odd point will lie on a line between the two adjacent even
points. The difference between this "prediction" and the actual odd
value replaces the odd element.
The update stage calculates the average of the odd and even
element pairs, although the method is indirect, since the
predict phase has over written the odd value.
Copyright and Use
You may use this source code without limitation and without
fee as long as you include:
This software was written and is copyrighted by Ian Kaplan, Bear
Products International, www.bearcave.com, 2001.
This software is provided "as is", without any warranty or
claim as to its usefulness. Anyone who uses this source code
uses it at their own risk. Nor is any support provided by
Ian Kaplan and Bear Products International.
Please send any bug fixes or suggested source changes to:
iank@bearcave.com
\author Ian Kaplan
*/
template
class line : public liftbase {
private:
/**
Calculate an extra "even" value for the line wavelet algorithm
at the end of the data series. Here we pretend that the
last two values in the data series are at the x-axis
coordinates 0 and 1, respectively. We then need to calculate
the y-axis value at the x-axis coordinate 2. This point lies
on a line running through the points at 0 and 1.
Given two points, x1, y1 and
x2, y2, where
x1 = 0
x2 = 1
calculate the point on the line at x3, y3,
where
x3 = 2
The "two-point equation" for a line given x1,
y1 and x2, y2 is
. y2 - y1
(y - y1) = -------- (x - x1)
. x2 - x1
Solving for y
. y2 - y1
y = -------- (x - x1) + y1
. x2 - x1
Since x1 = 0 and x2 = 1
. y2 - y1
y = -------- (x - 0) + y1
. 1 - 0
or
y = (y2 - y1)*x + y1
We're calculating the value at x3 = 2, so
y = 2*y2 - 2*y1 + y1
or
y = 2*y2 - y1
*/
double new_y( double y1, double y2)
{
double y = 2 * y2 - y1;
return y;
}
protected:
/**
Predict phase of line Lifting Scheme wavelet
The predict step attempts to "predict" the value of an
odd element from the even elements. The difference
between the prediction and the actual element is stored
as a wavelet coefficient.
The "predict" step takes place after the split step. The split
step will move the odd elements (bj) to the second half
of the array, leaving the even elements (ai) in the
first half
a0, a1, a1, a3, b0, b1, b2, b2,
The predict step of the line wavelet "predicts" that the
odd element will be on a line between two even elements.
bj+1,i = bj,i - (aj,i + aj,i+1)/2
Note that when we get to the end of the data series the odd
element is the last element in the data series (remember,
wavelet algorithms work on data series with 2n
elements). Here we "predict" that the odd element will be
on a line that runs through the last two even elements.
This can be calculated by assuming that the last two
even elements are located at x-axis coordinates 0 and 1,
respectively. The odd element will be at 2. The new_y()
function is called to do this simple calculation.
*/
void predict( T& vec, int N, transDirection direction )
{
int half = N >> 1;
double predictVal;
for (int i = 0; i < half; i++) {
int j = i + half;
if (i < half-1) {
predictVal = (vec[i] + vec[i+1])/2;
}
else if (N == 2) {
predictVal = vec[0];
}
else {
// calculate the last "odd" prediction
double n_plus1 = new_y( vec[i-1], vec[i] );
predictVal = (vec[i] + n_plus1)/2;
}
if (direction == forward) {
vec[j] = vec[j] - predictVal;
}
else if (direction == inverse) {
vec[j] = vec[j] + predictVal;
}
else {
printf("predictline::predict: bad direction value\n");
}
}
} // predict
/**
Update step of the linear interpolation wavelet
The predict phase works on the odd elements in the second
half of the array. The update phase works on the even
elements in the first half of the array. The update
phase attempts to preserve the average. After the update
phase is completed the average of the even elements should
be approximately the same as the average of the input data
set from the previous iteration. The result of the update
phase becomes the input for the next iteration.
In a Haar wavelet the average that replaces the even element is
calculated as the average of the even element and its associated
odd element (e.g., its odd neighbor before the split). This is
not possible in the line wavelet since the odd element has been
replaced by the difference between the odd element and the
mid-point of its two even neighbors. As a result, the odd element
cannot be recovered.
The value that is added to the even element to preserve the
average is calculated by the equation shown below. This
equation is given in Wim Sweldens' journal articles and his tutorial
(Building Your Own Wavelets at Home) and in
Ripples in Mathematics. A somewhat more complete derivation
of this equation is provided in Ripples in Mathematics
by A. Jensen and A. la Cour-Harbo, Springer, 2001.
The equation used to calculate the average is shown below for a
given iteratin i. Note that the predict phase has already
completed, so the odd values belong to iteration i+1.
eveni+1,j = eveni,j op (oddi+1,k-1 + oddi+1,k)/4
There is an edge problem here, when i = 0 and k = N/2 (e.g., there
is no k-1 element). We assume that the oddi+1,k-1 is
the same as oddk. So for the first element this
becomes
(2 * oddk)/4
or
oddk/2
*/
void update( T& vec, int N, transDirection direction )
{
int half = N >> 1;
for (int i = 0; i < half; i++) {
int j = i + half;
double val;
if (i == 0) {
val = vec[j]/2.0;
}
else {
val = (vec[j-1] + vec[j])/4.0;
}
if (direction == forward) {
vec[i] = vec[i] + val;
}
else if (direction == inverse) {
vec[i] = vec[i] - val;
}
else {
printf("update: bad direction value\n");
}
} // for
}
}; // line
#endif
packet/include/list.h 0100644 0000764 0000764 00000010243 07467260340 013543 0 ustar iank iank
#ifndef _LIST_H_
#define _LIST_H_
/** \file
The documentation in this file is formatted for doxygen
(see www.doxygen.org).
Copyright and Use
You may use this source code without limitation and without
fee as long as you include:
This software was written and is copyrighted by Ian Kaplan, Bear
Products International, www.bearcave.com, 2002.
This software is provided "as is", without any warranty or
claim as to its usefulness. Anyone who uses this source code
uses it at their own risk. Nor is any support provided by
Ian Kaplan and Bear Products International.
Please send any bug fixes or suggested source changes to:
iank@bearcave.com
@author Ian Kaplan
*/
/**
template class LIST
This is a generic list type. It should be instantiated
with a scalar type, like an integer or a pointer to
a larger type (e.g., a string or a structure). For example
LIST list;
LIST list;
The list "links" are allocated in a memory pool. This allocation
is handled by the local version of new in list_type. The memory
allocated for the list_type objects will be deallocated when
the memory pool is freed.
*/
template
class LIST
{
public:
/**
List backbone class
*/
class list_type {
public:
/** list data element */
T data;
/** pointer to the next element in the list */
list_type *next;
/** override the default new operator to allocate
list_type objects from the a memory pool */
void *operator new(size_t num_bytes)
{
block_pool mem_pool;
void *mem_addr = mem_pool.pool_alloc( num_bytes );
return mem_addr;
} // new
}; // class list_type
private:
/** list head */
list_type *list;
public:
/** define a handle type to abstract the list_type type */
typedef list_type *handle;
public:
/** class constructor */
LIST(void)
{
list = 0;
}
/** The copy constructor copies the list pointer */
LIST( const LIST &rhs )
{
list = rhs.list;
}
/** destructor does nothing */
~LIST() {}
/** deallocate the list */
void dealloc(void)
{
while ( remove() != 0 )
/* nada */;
} // dealloc
/** Add an element to the list */
void add( T data )
{
list_type *t;
t = new list_type();
t->data = data;
t->next = 0;
if (list == 0) {
list = t;
}
else {
t->next = list;
list = t;
}
} // add
/** reverse the list */
void reverse(void)
{
list_type *revlist = 0;
list_type *next;
for (list_type *t = list; t != 0; t = next) {
next = t->next;
t->next = revlist;
revlist = t;
}
list = revlist;
} // reverse
/** return the lenght of the list */
unsigned int length(void)
{
list_type *elem;
unsigned int cnt = 0;
for (elem = list; elem != 0; elem = elem->next)
cnt++;
return cnt;
} // lenght
/**
remove
Remove an element from the start of the list and return the first
element of the remaining list.
This function relies on the fact that the list elements were
allocated from a pool. The memory for these elements will be
recovered when the pool is deallocated.
*/
handle remove(void)
{
list_type *t;
if (list != 0) {
t = list;
list = t->next;
// no delete t;
}
return list;
} // remove
/** given a handle, return the associated data item */
T get_item( handle h)
{
return h->data;
} // get_item
/** get the first element from the list */
handle first(void)
{
return list;
} // first
/** iterator to get the next element */
handle next(handle h)
{
list_type *next = 0;
if (h != 0) {
next = h->next;
}
return next;
} // next
}; // template class LIST
#endif
packet/include/packcontainer.h 0100644 0000764 0000764 00000011300 10176075345 015404 0 ustar iank iank
#ifndef _PACKCONTAINER_H_
#define _PACKCONTAINER_H_
#include "packnode.h"
/**
The documentation in this file is formatted for doxygen
(see www.doxygen.org).
Copyright and Use
You may use this source code without limitation and without
fee as long as you include:
This software was written and is copyrighted by Ian Kaplan, Bear
Products International, www.bearcave.com, 2002.
This software is provided "as is", without any warranty or
claim as to its usefulness. Anyone who uses this source code
uses it at their own risk. Nor is any support provided by
Ian Kaplan and Bear Products International.
Please send any bug fixes or suggested source changes to:
iank@bearcave.com
@author Ian Kaplan
*/
/**
A container for use when calculating a packet wavelet tree.
By overriding the LHS and RHS versions of the [] operator,
packcontainer class allows a block of data to be treated as an array
while a wavelet transform step is being calculated. After the
wavelet transform step the data can be referenced as a left hand
side half (the result of the wavelet scaling function) and a right
hand side half (the result of the wavelet function). By allowing
the two halves of the array to be referenced, copying is avoided.
*/
class packcontainer {
private:
/** number of elements at this packet tree level */
size_t N;
/** left (low pass) half of the packcontainer data */
double* lhs;
/** right (high pass) half of the packnode data */
double* rhs;
private:
/** disallow the copy constructor */
packcontainer( const packcontainer &rhs ) {};
/** disallow default constructor */
packcontainer() {};
public:
/**
This version of teh packcontainer object is used as a container
in calculating the forward wavelet packet transform.
The packcontainer constructor is passed a pointer to a packnode
object. As the wavelet packet tree is being built this packnode
object contain the data from the previous level.
The length of the packnode object is N. The constructor will
dynamically allocate two vectors (lhs and rhs), each with N/2
elements. After the wavelet transform step has been calculated on
the packcontainer object, the lhs vector will contain the wavelet
scaling function, or low pass, result. The rhs vector will contain
the wavelet function, or high pass, result. These will be used
to construct two new packnode objects which will be children
of the object passed to the constructor.
*/
packcontainer( packnode* node )
{
assert( node != 0 );
N = node->length();
assert( N > 1 );
size_t half = N >> 1;
block_pool mem_pool;
size_t num_bytes = half * sizeof(double);
lhs = (double *)mem_pool.pool_alloc( num_bytes );
rhs = (double *)mem_pool.pool_alloc( num_bytes );
for (size_t i = 0; i < N; i++) {
(*this)[i] = (*node)[i];
}
} // packcontainer
/**
This version is used when calculating the inverse wavelet packet
transform. The constructor is passed the size of the container
(or at least the size that the container will be, when the left
and right hand size arrays are initialized). The lhs and rhs
arrays are then initialized and the inverse transform step
is calculated.
*/
packcontainer( size_t n )
{
N = n;
lhs = 0;
rhs = 0;
}
/**
This is a local implementation of new. When new
is applied to a packcontainer object, memory will
be allocated from a memory pool, rather than from
the system memory pool.
*/
void *operator new( unsigned int num_bytes )
{
block_pool mem_pool;
void *mem_addr = mem_pool.pool_alloc( num_bytes );
return mem_addr;
} // new
/** LHS [] operator */
double &operator[]( const size_t i )
{
assert( i < N );
size_t half = N >> 1;
if (i < half)
return lhs[i];
else {
return rhs[i-half];
}
}
/** RHS [] operator */
double operator[]( const size_t i ) const
{
assert( i < N );
size_t half = N >> 1;
if (i < half)
return lhs[i];
else {
return rhs[i-half];
}
}
/** return the left hand size array */
double* lhsData() { return lhs; }
/** return the right hand size array */
double* rhsData() { return rhs; }
/** set the left hand size array */
void lhsData(double* l) { lhs = l; }
/** set the right hand size array */
void rhsData(double* r) { rhs = r; }
/** return the length of the data in the packet
container. Note that this length is
the length of the rhs plus the length of
the lhs arrays.
*/
size_t length() { return N; }
}; // packcontainer
#endif
packet/include/packdata.h 0100644 0000764 0000764 00000006126 07525053412 014340 0 ustar iank iank
#ifndef _PACKDATA_H_
#define _PACKDATA_H_
/** \file
The documentation in this file is formatted for doxygen
(see www.doxygen.org).
Copyright and Use
You may use this source code without limitation and without
fee as long as you include:
This software was written and is copyrighted by Ian Kaplan, Bear
Products International, www.bearcave.com, 2002.
This software is provided "as is", without any warranty or
claim as to its usefulness. Anyone who uses this source code
uses it at their own risk. Nor is any support provided by
Ian Kaplan and Bear Products International.
Please send any bug fixes or suggested source changes to:
iank@bearcave.com
@author Ian Kaplan
*/
#include "blockpool.h"
/**
The packdata class is a container for the core information
for a wavelet packet transform node. This class serves
as the base class for the packnode class. It is also
the element used in the list that results from the wavelet
packet transform. This list represents the data in a minimal
for, relative to the wavelet packet transform cost function.
\author Ian Kaplan
*/
template
class packdata {
private:
/** Disallow the copy constructor */
packdata( const packdata &rhs ) {}
public:
typedef enum { BadNodeKind,
OriginalData,
LowPass,
HighPass } transformKind;
protected:
/** Kind of data: original data, result of the
low pass filer (wavelet scaling function), result
of the high pass filter (wavelet function) */
transformKind kind;
/** number of elements at this packet tree level */
size_t N;
/** Wavelet packet data */
T* data;
/** default constructor */
packdata() {}
public:
/**
\arg vec a pointer to an array of type T elements
\arg n the size of the array of type T elements
\arg k the kind of data (e.g., original data, low pass result,
high pass result.
*/
packdata( T *vec, const size_t n, const transformKind k )
{
data = vec;
N = n;
kind = k;
}
/** define a virtual destructure which in the base class */
virtual ~packdata() {}
/**
Overload the standard new operator and allocate
memory from the memory pool.
*/
void *operator new( unsigned int num_bytes )
{
block_pool mem_pool;
void *mem_addr = mem_pool.pool_alloc( num_bytes );
return mem_addr;
} // new
/** print the data */
void pr() const
{
for (int i = 0; i < N; i++) {
printf("%7.4f ", data[i] );
}
printf("\n");
} // pr
/** get a pointer to the data */
const T* getData() { return data; }
/**
Return the "kind" of data (e.g., original data,
low pass result, high pass result).
*/
const transformKind getKind() { return kind; }
/** return the length of the array in the packdata object */
size_t length() { return N; }
}; // packdata
#endif
packet/include/packdata_list.h 0100644 0000764 0000764 00000003564 07525053460 015401 0 ustar iank iank
#ifndef _PACKDATA_LIST_H_
#define _PACKDATA_LIST_H_
#include "fifo_list.h"
#include "packdata.h"
/** \file
The documentation in this file is formatted for doxygen
(see www.doxygen.org).
Copyright and Use
You may use this source code without limitation and without
fee as long as you include:
This software was written and is copyrighted by Ian Kaplan, Bear
Products International, www.bearcave.com, 2002.
This software is provided "as is", without any warranty or
claim as to its usefulness. Anyone who uses this source code
uses it at their own risk. Nor is any support provided by
Ian Kaplan and Bear Products International.
Please send any bug fixes or suggested source changes to:
iank@bearcave.com
\author Ian Kaplan
*/
/**
The packdata_list class is based on a subclass built from a
FIFO_LIST instance of the FIFO_LIST template.
The packdata_list is constructed by traversing the wavelet packet
tree depth first, left to right. New elements are added to the end
of the list, ordering the list from front to back. The resulting
list reflects the ordering of the wavelet packet tree nodes the make
up the best basis.
This class extends the FIFL_LIST instance by adding a print function.
\author Ian Kaplan
*/
template
class packdata_list : public FIFO_LIST *>
{
public:
/**
Print the packet data list. Each list element (which is a wavelet
packet data set) is printed one element per line.
*/
void pr()
{
handle h;
packdata *elem;
for (h = first(); h != 0; h = next( h )) {
elem = get_item( h );
elem->pr();
} // for h
printf("\n");
} // pr
}; // packdata_list
#endif
packet/include/packfreq.h 0100644 0000764 0000764 00000005166 07525053515 014373 0 ustar iank iank
#ifndef _PACKFREQ_H_
#define _PACKFREQ_H_
#include "packnode.h"
#include "packcontainer.h"
#include "packtree_base.h"
#include "liftbase.h"
#include "grow_array.h"
/** \file
The documentation in this file is formatted for doxygen
(see www.doxygen.org).
Copyright and Use
You may use this source code without limitation and without
fee as long as you include:
This software was written and is copyrighted by Ian Kaplan, Bear
Products International, www.bearcave.com, 2002.
This software is provided "as is", without any warranty or
claim as to its usefulness. Anyone who uses this source code
uses it at their own risk. Nor is any support provided by
Ian Kaplan and Bear Products International.
Please send any bug fixes or suggested source changes to:
iank@bearcave.com
\author Ian Kaplan
*/
/**
Build a wavelet packet tree for frequency analysis.
Wavelet frequency analysis uses a modified wavelet packet tree.
Horizontal slices through the modified wavelet packet tree (the so
called "level basis") are ordered in increasing frequency regions.
This class and the standard wavelet packet tree class are derived
from the same base class. The standard wavelet class includes
functions to calculate a minimal data representation relative
to a cost function.
The constructor for this class is passed two arguments:
-
A vector of doubles containing the data set (e.g., the signal).
The length of the vector must be a power of two.
-
A pointer to a wavelet lifting scheme class that will be
used in calculating the wavelet transform step.
If the vector passed to the constructor contains N double values, the
result of the constructor will be a wavelet packet tree with
log2(N) levels.
\author Ian Kaplan
*/
class packfreq : public packtree_base {
private:
/** Level basis matrix */
GrowableArray *> mat;
void findLevel( packnode* top,
size_t cur_level,
const size_t level );
protected:
/** disallow the copy constructor */
packfreq( const packfreq &rhs ) {};
/** disallow the default constructor */
packfreq() {};
public:
packfreq( const double *vec,
const size_t n,
liftbase *w );
/** destructor does nothing */
~packfreq() {}
void getLevel( const size_t level );
void plotMat(const size_t N);
void prMat();
}; // packfreq
#endif
packet/include/packnode.h 0100644 0000764 0000764 00000007444 07521557425 014371 0 ustar iank iank
#ifndef _PACKNODE_H_
#define _PACKNODE_H_
#include "packdata.h"
/** \file
The documentation in this file is formatted for doxygen
(see www.doxygen.org).
Copyright and Use
You may use this source code without limitation and without
fee as long as you include:
This software was written and is copyrighted by Ian Kaplan, Bear
Products International, www.bearcave.com, 2002.
This software is provided "as is", without any warranty or
claim as to its usefulness. Anyone who uses this source code
uses it at their own risk. Nor is any support provided by
Ian Kaplan and Bear Products International.
Please send any bug fixes or suggested source changes to:
iank@bearcave.com
@author Ian Kaplan
*/
/**
A wavelet packet tree node
A description of the wavelet packet algorithm can be found
in Chapter 8 of Ripples in Mathematics by Jensen and
la Cour-Harbo.
For a data set consisting of N data elements, the wavelet packet
algorithm will build a binary tree with log2(N) levels.
Since it is a binary tree, the number of nodes doubles at each
level. At the same time, the amount of data that is stored in
each node is halved.
A node will have n data elements. The two children will
each have n/2 elements, stored in packnode children.
A pointer to the data vector for the node and the length of the
data is set by the class constructor.
Once the wavelet packet tree is built a cost value is assigned to
each node in the tree.
*/
template
class packnode : public packdata {
private:
/** left child (with N/2) data elements */
packnode* leftChild;
/** right child (with N/2) data elements */
packnode* rightChild;
/** cost value for this level */
T costVal;
/** chosen == true: node is part of the best basis of the
wavelet transform, otherwise, false. */
bool chosen;
private:
/** disallow the copy constructor */
packnode( const packnode &rhs ) {};
/** disallow default constructor */
packnode() {};
public:
/**
Packnode constructor
*/
packnode( T *vec,
const size_t n,
const transformKind k ) : packdata(vec, n, k)
{
leftChild = 0;
rightChild = 0;
costVal = 0.0;
chosen = false;
}
/** LHS [] operator */
T &operator[]( const size_t i )
{
assert( i < N );
return data[i];
}
/** RHS [] operator */
T operator[]( const size_t i ) const
{
assert( i < N );
return data[i];
}
/** print the cost value */
void prCost() const
{
printf("%7.4f\n", costVal );
}
/** if the node is selected as part of the best basis
function, print it.
*/
void prBestBasis() const
{
for (int i = 0; i < N; i++) {
printf("%7.4f ", data[i] );
}
if (chosen) {
printf(" *");
}
printf("\n");
} // prBestBasis
/** set the left child pointer */
void lhsChild( packnode *l ) { leftChild = l; }
/** get the left child pointer */
packnode *lhsChild(void) { return leftChild; }
/** set the right child pointer */
void rhsChild( packnode *r ) { rightChild = r; }
/** get the right child pointer */
packnode *rhsChild(void) { return rightChild; }
/** set the cost value for the node */
void cost( T val ) { costVal = val; }
/** get the cost value for the node */
T cost(void) { return costVal; }
/** the "chosen" flag marks a node for inclusion in
the best basis set.
*/
void mark( bool b ) { chosen = b; }
/** return the value of the "mark" boolean flag */
bool mark() { return chosen; }
}; // packnode
#endif
packet/include/packtree.h 0100644 0000764 0000764 00000005011 07515163704 014363 0 ustar iank iank
#ifndef _PACKTREE_H_
#define _PACKTREE_H_
#include "packtree_base.h"
#include "packdata_list.h"
#include "packcontainer.h"
#include "liftbase.h"
/** \file
The documentation in this file is formatted for doxygen
(see www.doxygen.org).
Copyright and Use
You may use this source code without limitation and without
fee as long as you include:
This software was written and is copyrighted by Ian Kaplan, Bear
Products International, www.bearcave.com, 2002.
This software is provided "as is", without any warranty or
claim as to its usefulness. Anyone who uses this source code
uses it at their own risk. Nor is any support provided by
Ian Kaplan and Bear Products International.
Please send any bug fixes or suggested source changes to:
iank@bearcave.com
@author Ian Kaplan
*/
/**
The packtree object constructs a wavelet packet tree
The constructor is passed a vector of doubles, the length of the
vector (which must be a power of two) and a pointer to a wavelet
lifting scheme class that will be used in calculating the wavelet
transform step. If the vector passed to the constructor contains N
double vaues, the result of the constructor will be a wavelet packet
tree with log2(N) levels.
*/
class packtree : public packtree_base {
private:
/** Found original data marked as part of the best basis.
This means that the best basis function failed (or
that the original data is the most compact representation
relative to the cost function used). */
bool foundOriginalData;
/** found a best basis value in the wavelet packet tree */
bool foundBestBasisVal;
private:
/** disallow the copy constructor */
packtree( const packtree &rhs ) {};
/** disallow the default constructor */
packtree() {};
double bestBasisWalk( packnode *root );
void buildBestBasisList( packnode *root,
packdata_list &list );
void checkBestBasis( packnode *root );
void cleanTree(packnode *root, bool removeMark );
public:
packtree( const double *vec,
const size_t n,
liftbase *w );
/** destructor does nothing */
~packtree() {}
void prCost();
void prBestBasis();
void bestBasis();
bool bestBasisOK();
packdata_list getBestBasisList();
}; // packtree
#endif
packet/include/packtree_base.h 0100644 0000764 0000764 00000003175 10176075625 015367 0 ustar iank iank
#ifndef _PACKTREE_BASE_H_
#define _PACKTREE_BASE_H_
/** \file
The documentation in this file is formatted for doxygen
(see www.doxygen.org).
Copyright and Use
You may use this source code without limitation and without
fee as long as you include:
This software was written and is copyrighted by Ian Kaplan, Bear
Products International, www.bearcave.com, 2002.
This software is provided "as is", without any warranty or
claim as to its usefulness. Anyone who uses this source code
uses it at their own risk. Nor is any support provided by
Ian Kaplan and Bear Products International.
Please send any bug fixes or suggested source changes to:
iank@bearcave.com
\author Ian Kaplan
*/
#include "packnode.h"
#include "liftbase.h"
/**
Base class for wavelet packet trees. Subclasses for this
class include classes to built wavelet packet trees for
best basis calculation and wavelet packet trees for
frequency analysis.
\author Ian Kaplan
*/
class packtree_base {
protected:
/** root of the wavelet packet tree */
packnode *root;
/** wavelet packet transform object */
liftbase *waveObj;
typedef enum { BadPrintKind,
printData,
printCost,
printBestBasis } printKind;
void breadthFirstPrint(printKind kind);
void newLevel( packnode* top, bool freqCalc, bool reverse );
public:
void pr();
/** get the root of the wavelet packet tree */
packnode *getRoot() { return root; }
}; // packtree_base
#endif
packet/include/queue.h 0100644 0000764 0000764 00000006076 07525053631 013723 0 ustar iank iank
#ifndef _QUEUE_H_
#define _QUEUE_H_
#include "fifo_list.h"
#include "packnode.h"
/** \file
This file defines a queue (FIFO) object for use in
breadth first tree traversal of a wavelet packet tree.
See chapter 8 of Ripples in Mathematics by
Jensen and la Cour-Harbo for a description of
the wavelet packet algorithm.
The documentation in this file is formatted for doxygen
(see www.doxygen.org).
Copyright and Use
You may use this source code without limitation and without
fee as long as you include:
This software was written and is copyrighted by Ian Kaplan, Bear
Products International, www.bearcave.com, 2002.
This software is provided "as is", without any warranty or
claim as to its usefulness. Anyone who uses this source code
uses it at their own risk. Nor is any support provided by
Ian Kaplan and Bear Products International.
Please send any bug fixes or suggested source changes to:
iank@bearcave.com
@author Ian Kaplan
*/
/**
A queue element. This class stores a pointer to a wavelet
packet tree node (a packnode) and the indentation (in spaces)
to be used when printing the tree.
\author Ian Kaplan
*/
template
class queueElem
{
private:
/** disallow the default constructor */
queueElem() {}
public:
packnode *node;
size_t indent;
/** Initialize the queue element */
queueElem( packnode *n, size_t i )
{
node = n;
indent = i;
}
/** allocate queueElem objects from a memory pool */
void *operator new(size_t num_bytes)
{
block_pool mem_pool;
void *mem_addr = mem_pool.pool_alloc( num_bytes );
return mem_addr;
} // new
}; // queueElem
/**
The queue class extends the FIFO_LIST template (which
in this case is instantiated with queueElem *. This
class is designed to support breadth first printing
of a wavelet packet tree.
\author Ian Kaplan
*/
template
class queue : protected FIFO_LIST *>
{
public:
/** Get the first element in the queue */
queueElem *queueStart()
{
handle h = first();
queueElem *elem = get_item( h );
return elem;
} // queueStart
/**
Remove the element at the start of the list. This function does
not actually call delete to recover the object. It relies on the
fact that the FIFO_LIST template uses pool allocation and the
memory will be recovered when the pool is deallocated.
*/
void deleteStart()
{
handle h = first();
if (h != 0) {
queueElem *elem = get_item( h );
remove();
// no delete elem;
}
} // deleteStart
/** Add an element to the queue */
void addQueue(packnode *node, size_t indent )
{
queueElem *elem = new queueElem(node, indent);
add( elem );
} // addQueue
/** return true if the queue is empty */
bool queueEmpty() { return (first() == 0); }
}; // queue
#endif
packet/include/packcontainer.h~ 0100644 0000764 0000764 00000011570 07515162100 015600 0 ustar iank iank
#ifndef _PACKCONTAINER_H_
#define _PACKCONTAINER_H_
#include "packnode.h"
/**
The documentation in this file is formatted for doxygen
(see www.doxygen.org).
Copyright and Use
You may use this source code without limitation and without
fee as long as you include:
This software was written and is copyrighted by Ian Kaplan, Bear
Products International, www.bearcave.com, 2002.
This software is provided "as is", without any warranty or
claim as to its usefulness. Anyone who uses this source code
uses it at their own risk. Nor is any support provided by
Ian Kaplan and Bear Products International.
Please send any bug fixes or suggested source changes to:
iank@bearcave.com
@author Ian Kaplan
*/
/**
A container for use when calculating a packet wavelet tree.
By overriding the LHS and RHS versions of the [] operator,
packcontainer class allows a block of data to be treated as an array
while a wavelet transform step is being calculated. After the
wavelet transform step the data can be referenced as a left hand
side half (the result of the wavelet scaling function) and a right
hand side half (the result of the wavelet function). By allowing
the two halves of the array to be referenced, copying is avoided.
*/
class packcontainer {
private:
/** number of elements at this packet tree level */
size_t N;
/** left (low pass) half of the packcontainer data */
double* lhs;
/** right (high pass) half of the packnode data */
double* rhs;
private:
/** disallow the copy constructor */
packcontainer( const packcontainer &rhs ) {};
/** disallow default constructor */
packcontainer() {};
public:
/**
This version of teh packcontainer object is used as a container
in calculating the forward wavelet packet transform.
The packcontainer constructor is passed a pointer to a packnode
object. As the wavelet packet tree is being built this packnode
object contain the data from the previous level.
The length of the packnode object is N. The constructor will
dynamically allocate two vectors (lhs and rhs), each with N/2
elements. After the wavelet transform step has been calculated on
the packcontainer object, the lhs vector will contain the wavelet
scaling function, or low pass, result. The rhs vector will contain
the wavelet function, or high pass, result. These will be used
to construct two new packnode objects which will be children
of the object passed to the constructor.
*/
packcontainer( packnode* node )
{
assert( node != 0 );
N = node->length();
assert( N > 1 );
size_t half = N >> 1;
block_pool mem_pool;
size_t num_bytes = half * sizeof(double);
lhs = (double *)mem_pool.pool_alloc( num_bytes );
rhs = (double *)mem_pool.pool_alloc( num_bytes );
for (size_t i = 0; i < N; i++) {
(*this)[i] = (*node)[i];
}
} // packcontainer
/**
This version is used when calculating the inverse wavelet packet
transform. The constructor is passed the size of the container
(or at least the size that the container will be, when the left
and right hand size arrays are initialized). The lhs and rhs
arrays are then initialized and the inverse transform step
is calculated.
*/
packcontainer( size_t n )
{
N = n;
lhs = 0;
rhs = 0;
}
/**
This is a local implementation of new. When new
is applied to a packcontainer object, memory will
be allocated from a memory pool, rather than from
the system memory pool.
*/
void *operator new( unsigned int num_bytes )
{
block_pool mem_pool;
void *mem_addr = mem_pool.pool_alloc( num_bytes );
return mem_addr;
} // new
/** LHS [] operator */
double &operator[]( const size_t i )
{
assert( i < N );
size_t half = N >> 1;
if (i < half)
return lhs[i];
else {
return rhs[i-half];
}
}
/** RHS [] operator */
double operator[]( const size_t i ) const
{
assert( i < N );
size_t half = N >> 1;
if (i < half)
return lhs[i];
else {
return rhs[i-half];
}
}
/** return the left hand size array */
double* lhsData() { return lhs; }
/** return the right hand size array */
double* rhsData() { return rhs; }
/** set the left hand size array */
void lhsData(double* l) { lhs = l; }
/** set the right hand size array */
void rhsData(double* r) { rhs = r; }
/** return the length of the data in the packet
container. Note that this length is
the length of the rhs plus the length of
the lhs arrays.
*/
size_t length() { return N; }
}; // packcontainer
#endif
packet/include/packtree_base.h~ 0100644 0000764 0000764 00000003310 07525053573 015556 0 ustar iank iank
#ifndef _PACKTREE_BASE_H_
#define _PACKTREE_BASE_H_
/** \file
The documentation in this file is formatted for doxygen
(see www.doxygen.org).
Copyright and Use
You may use this source code without limitation and without
fee as long as you include:
This software was written and is copyrighted by Ian Kaplan, Bear
Products International, www.bearcave.com, 2002.
This software is provided "as is", without any warranty or
claim as to its usefulness. Anyone who uses this source code
uses it at their own risk. Nor is any support provided by
Ian Kaplan and Bear Products International.
Please send any bug fixes or suggested source changes to:
iank@bearcave.com
\author Ian Kaplan
*/
#include "packnode.h"
#include "liftbase.h"
/**
Base class for wavelet packet trees. Subclasses for this
class include classes to built wavelet packet trees for
best basis calculation and wavelet packet trees for
frequency analysis.
\author Ian Kaplan
*/
class packtree_base {
protected:
/** root of the wavelet packet tree */
packnode *root;
/** wavelet packet transform object */
liftbase *waveObj;
typedef enum { BadPrintKind,
printData,
printCost,
printBestBasis } printKind;
void breadthFirstPrint(printKind kind);
void newLevel( packnode* top, bool freqCalc, bool reverse );
public:
void pr();
/** get the root of the wavelet packet tree */
packnode *getRoot() { return root; }
}; // packtree_base
#endif
packet/lossless/ 0040755 0000764 0000764 00000000000 07525263472 012652 5 ustar iank iank packet/lossless/include/ 0040755 0000764 0000764 00000000000 07525060126 014264 5 ustar iank iank packet/lossless/include/costbase_int.h 0100644 0000764 0000764 00000005141 07523327722 017116 0 ustar iank iank
#ifndef _COSTBASE_INT_H_
#define _COSTBASE_INT_H_
#include "packnode.h"
/** \file
The documentation in this file is formatted for doxygen
(see www.doxygen.org).
Copyright and Use
You may use this source code without limitation and without
fee as long as you include:
This software was written and is copyrighted by Ian Kaplan, Bear
Products International, www.bearcave.com, 2002.
This software is provided "as is", without any warranty or
claim as to its usefulness. Anyone who uses this source code
uses it at their own risk. Nor is any support provided by
Ian Kaplan and Bear Products International.
Please send any bug fixes or suggested source changes to:
iank@bearcave.com
@author Ian Kaplan
*/
/**
Base class for objects that define integer costs functions
for an integer version of the wavelet packet transform.
The costbase base class provides a constructor that is passed the
root of an integer wavelet packet tree (which is built by
packtree_int.cpp). The integer wavelet packet tree is constructed
from packnode objects (which are a subclass of the packdata).
The cost function calculation invoked by the constructor traverses
the wavelet packet tree (top down) and calls costCalc on each node
in the tree. The cost function result is stored in the node. Note
that the cost function also calculates a cost value for the original
data, since in theory the original data may represent the minimal
representation for the data in terms of the cost function.
The pure virtual function costCalc, which calculates the cost
function, must be defined by the subclass.
A description of the cost functions associated with the wavelet
packet transform can be found in Chapter 8 of Ripples in
Mathematics by Jense and la Cour-Harbo.
*/
class costbase_int
{
private:
/** disallow the copy constructor */
costbase_int( const costbase_int &rhs ) {}
protected:
/**
Recursively traverse the wavelet packet tree and calculate the cost
function.
*/
void traverse( packnode *node )
{
if (node != 0) {
int cost = costCalc( node );
node->cost( cost );
traverse( node->lhsChild() );
traverse( node->rhsChild() );
}
} // traverse
/** Cost function to be defined by the subclass */
virtual int costCalc(packnode *node) = 0;
public:
/** The default constructor does nothing */
costbase_int() {}
}; // costbase_int
#endif
packet/lossless/include/costwidth.h 0100644 0000764 0000764 00000002354 07523334101 016441 0 ustar iank iank
#ifndef _COSTWIDTH_H_
#define _COSTWIDTH_H_
/** \file
The documentation in this file is formatted for doxygen
(see www.doxygen.org).
Copyright and Use
You may use this source code without limitation and without
fee as long as you include:
This software was written and is copyrighted by Ian Kaplan, Bear
Products International, www.bearcave.com, 2002.
This software is provided "as is", without any warranty or
claim as to its usefulness. Anyone who uses this source code
uses it at their own risk. Nor is any support provided by
Ian Kaplan and Bear Products International.
Please send any bug fixes or suggested source changes to:
iank@bearcave.com
@author Ian Kaplan
*/
#include "costbase_int.h"
/**
For a wavelet packet tree node, consisting of n integer
values, calculate the minimal number of bits needed to represent
the values. This cost function is used for lossless wavelet
packet compression.
*/
class costwidth : public costbase_int
{
protected:
int costCalc( packnode *root );
public:
costwidth( packnode *root ) { traverse( root ); }
};
#endif
packet/lossless/include/delta.h 0100644 0000764 0000764 00000003544 07525051002 015522 0 ustar iank iank
/**
The delta transform.
A data set of N elements is replaced by the differences. In this
simple algorithm
si = si - si-1.
For an array indexed from 0, i = 1 to N-1. The element at s0
is the reference element, which is unchanged. The next element
at s1 is replaced by the difference between that element
and s0.
The algorithm implemented here is an in-place algorithm which
replaces the original data in the forward transform and
reconstructs the data from in-place in the inverse transform.
The class is implemented as a template class. Obviously the
type used to instantiate the template must support signed
arithmetic.
The class was written as a comparision baseline for the wavelet
compresson algorithms. The "delta" algorithm is very simple and has
a time complexity of N, where as the wavelet algorithm is
Nlog2N. Obviously, if the wavelet algorithm does
not do better than this simple algorithm for the data set in
question, the wavelet is a poor choice (or, perhaps, the wavelet
function in the predict step is poorly chosen for the data set).
*/
template
class delta_trans
{
public:
/**
Convert the value in the array into an initial reference
value and a set of delta values.
*/
void forward( T *vec, size_t len )
{
if (vec != 0 && len > 0) {
// reference value
T next;
T refVal = vec[0];
for (size_t i = 1; i < len; i++) {
next = vec[i];
vec[i] = vec[i] - refVal;
refVal = next;
}
}
} // forward
void inverse( T *vec, size_t len )
{
if (vec != 0 && len > 0) {
for (size_t i = 1; i < len; i++) {
vec[i] = vec[i] + vec[i-1];
}
}
} // inverse
}; // delta_trans
packet/lossless/include/haar_int.h 0100644 0000764 0000764 00000011075 07474537664 016246 0 ustar iank iank #ifndef _HAAR_INT_H_
#define _HAAR_INT_H_
/** \file
The documentation in this file is formatted for doxygen
(see www.doxygen.org).
Copyright and Use
You may use this source code without limitation and without
fee as long as you include:
This software was written and is copyrighted by Ian Kaplan, Bear
Products International, www.bearcave.com, 2002.
This software is provided "as is", without any warranty or
claim as to its usefulness. Anyone who uses this source code
uses it at their own risk. Nor is any support provided by
Ian Kaplan and Bear Products International.
Please send any bug fixes or suggested source changes to:
iank@bearcave.com
@author Ian Kaplan
*/
#include "stdio.h"
#include "liftbase.h"
/**
A lifting scheme version of the Haar integer to integer transform.
The standard wavelet transform creates real wavelet coefficients,
even if the input data consists of integers. This is a problem in
lossless compression (e.g., lossless image compression) and in other
compression related algorithm.
This verson of the Haar wavelet transform takes an data set and
creates an integer result. In the case of the Lifting Scheme
version of the Haar transform, the code is the same as the
real version, except that integers are used.
This algorithm is sometimes called the S-transform in the
image compression world.
References
-
Wavelet Transforms that Map Integers to Integers by
A.R. Calderbank, ingrid daubechies, wim weldens and Boon-Lock Yeo,
August 1996
This is the central reference that was used to develop this code.
Parts 1 and 2 of this paper are for the mathematicially
sophisticated (which is to say, they are not light reading).
However, for the implementer, part 3 and part 4 of this paper
provide excellent coverage of perfectly invertable wavelet
transforms that map integers to integers. In fact, part 3 of this
paper is worth reading in general for its discussion of the wavelet
Lifting Scheme.
-
Ripples in Mathematics: the Discrete Wavelet Transform
by Arne Jense and Anders la Cour-Harbo, Springer, 2001
This book is a good reference for the Lifting Scheme and the
wavelet transform in general.
*/
class haar_int : public liftbase
{
public:
/** the constructor does nothing */
haar_int() {}
/** the destructor does nothing */
~haar_int() {}
/** declare but do not define the copy constructor */
haar_int( const haar_int &rhs );
protected:
/**
Haar wavelet lifting scheme predict step.
The predict step "predicts" that an odd value will be
equal to the even value. The difference between the
actual value of the odd element and the even element
are stored in the upper half of the array.
The predict step is sometime referred to as the high
pass filter or the wavelet function.
The integer wavelet transform predict step is the same
as the standard (real) version of the lifting scheme
Haar transform.
*/
void predict( int *& vec, int N, transDirection direction )
{
int half = N >> 1;
for (int i = 0; i < half; i++) {
int predictVal = vec[i];
int j = i + half;
if (direction == forward) {
vec[j] = vec[j] - predictVal;
}
else if (direction == inverse) {
vec[j] = vec[j] + predictVal;
}
else {
printf("haar_int::predict: bad direction value\n");
}
}
} // predict
/**
Update step of the integer to integer wavelet transform.
In the Haar transform the update step calculates the low pass
filter (or average). For a detailed discussion of this
algorithm, see
Basic Lifting Scheme Wavelets.
*/
void update( int *& vec, int N, transDirection direction )
{
int half = N >> 1;
for (int i = 0; i < half; i++) {
int j = i + half;
// updateVal = floor( vec[j] / 2.0 )
int updateVal = vec[j] >> 1;
if (direction == forward) {
vec[i] = vec[i] + updateVal;
}
else if (direction == inverse) {
vec[i] = vec[i] - updateVal;
}
else {
printf("update_int: bad direction value\n");
}
}
} // update
}; // haar_int
#endif
packet/lossless/include/invpacktree_int.h 0100644 0000764 0000764 00000006351 07523334250 017624 0 ustar iank iank
#ifndef _INVPACKTREE_INT_H_
#define _INVPACKTREE_INt_H_
/** \file
The documentation in this file is formatted for doxygen
(see www.doxygen.org).
Copyright and Use
You may use this source code without limitation and without
fee as long as you include:
This software was written and is copyrighted by Ian Kaplan, Bear
Products International, www.bearcave.com, 2002.
This software is provided "as is", without any warranty or
claim as to its usefulness. Anyone who uses this source code
uses it at their own risk. Nor is any support provided by
Ian Kaplan and Bear Products International.
Please send any bug fixes or suggested source changes to:
iank@bearcave.com
@author Ian Kaplan
*/
#include "liftbase.h"
#include "list.h"
#include "packcontainer_int.h"
#include "packdata.h"
#include "packdata_list.h"
/**
Inverse wavelet packet transform
The invpacktree_int constructor is passed a packdata_list object and a
wavelet transform object. It calculates the inverse wavelet
packet transform, using the data in the packdata_list object
and the inverse wavelet transform step function of the wavelet
transform object. The best basis data is destroyed in calculating
the inverse transform.
The packdata_list object contains the "best basis" result from a
wavelet packet transform. The wavelet packet transform should have
been calculated with the same wavelet transform as the object passed
to this constructor.
After the constructor completes, the data result can be
obtained by calling the getData() function.
The wavelet transforms used by this object are all derived
from the liftbase class and are "lifting scheme" wavelet
transforms.
I have found the wavelet literature difficult, in general. When it
comes to the wavelet packet transform I have found most of it
impossible to understand. Impossible, that it until I got the book
Ripples in Mathematics by Jensen and la Cour-Harbo, Springer
Verlag, 2001. The wavelet packet transform, for which this class
is the inverse, is heavily based on Chapter 8 of Ripples in
Mathematics.
I have found very little material on the actual implementation of
the inverse wavelet packet transform. This algorithm is my own
design.
*/
class invpacktree_int {
private:
/** wavelet transform object */
liftbase *waveObj;
/** disallow the copy constructor */
invpacktree_int( const invpacktree_int &rhs ) {}
/** inverse wavelet packet transform calculation stack */
LIST stack;
/** pointer to the inverse transform result */
const int *data;
/** length of data */
size_t N;
private:
void new_level( packdata *elem );
void add_elem( packdata *elem );
void reduce();
public:
invpacktree_int( packdata_list &list,
liftbase *w );
/** The destructor does nothing */
~invpacktree_int() {}
/** Get the result of the inverse packet transform */
const int *getData() { return data; }
void pr();
}; // invpacktree_int
#endif
packet/lossless/include/line_int.h 0100644 0000764 0000764 00000021345 07523333620 016240 0 ustar iank iank
#ifndef _LINE_INT_H_
#define _LINE_INT_H_
/** \file
The documentation in this file is formatted for doxygen
(see www.doxygen.org).
Copyright and Use
You may use this source code without limitation and without
fee as long as you include:
This software was written and is copyrighted by Ian Kaplan, Bear
Products International, www.bearcave.com, 2002.
This software is provided "as is", without any warranty or
claim as to its usefulness. Anyone who uses this source code
uses it at their own risk. Nor is any support provided by
Ian Kaplan and Bear Products International.
Please send any bug fixes or suggested source changes to:
iank@bearcave.com
@author Ian Kaplan
*/
#include "stdio.h"
#include "liftbase.h"
/**
An integer version of the linear interpolation wavelet
The linear interpolation wavelet uses a predict phase
that "predicts" that an odd element in the data set
will line on a line between its two even neighbors.
This is an integer version of the linear interpolation wavelet. It
is interesting to note that unlike the S transform (the integer
version of the Haar wavelet) or the TS transform (an integer version
of the CDF(3,1) transform) this algorithm does not preserve
the mean. That is, when the transform is calculated, the first
element of the result array will not be the mean.
*/
template
class line_int : public liftbase
{
public:
/** the constructor does nothing */
line_int() {}
/** the destructor does nothing */
~line_int() {}
/** declare, but do not define the copy constructor */
line_int( const line_int &rhs );
private:
/**
Given y1 at x-coordinate 0 and y2 at x-coordinate
1, calculate y, at x-coordinate 2.
*/
int new_n_plus1( int y1, int y2)
{
int y = 2 * y2 - y1;
return y;
}
/**
Given a point y1 at x-coordinate 0 and y2 at x-coordinate 1,
calculate y at x-coordinate -1.
*/
int new_n_minus1( int y1, int y2)
{
int y = 2 * y1 - y2;
return y;
}
protected:
/**
Predict phase of Lifting Scheme linear interpolation wavelet
The predict step attempts to "predict" the value of an
odd element from the even elements. The difference
between the prediction and the actual element is stored
as a wavelet coefficient.
The "predict" step takes place after the split step. The split
step will move the odd elements (bj) to the second half
of the array, leaving the even elements (ai) in the
first half
a0, a1, a1, a3, b0, b1, b2, b2,
The predict step of the line wavelet "predicts" that the
odd element will be on a line between two even elements.
bj+1,i = bj,i - (aj,i + aj,i+1)/2
Note that when we get to the end of the data series the odd
element is the last element in the data series (remember, wavelet
algorithms work on data series with 2n elements). Here
we "predict" that the odd element will be on a line that runs
through the last two even elements. This can be calculated by
assuming that the last two even elements are located at x-axis
coordinates 0 and 1, respectively. The odd element will be at 2.
The new_n_plus1() function is called to do this simple
calculation.
Note that in the case where (N == 2), the algorithm becomes
the same as the Haar wavelet. We "predict" that the odd value
vec[1] will be the same as the even value, vec[0].
*/
void predict( T & vec, int N, transDirection direction )
{
int half = N >> 1;
int predictVal;
for (int i = 0; i < half; i++) {
int j = i + half;
if (i < half-1) {
predictVal = (int)((((float)vec[i] + (float)vec[i+1])/2.0) + 0.5);
}
else if (N == 2) {
predictVal = vec[0];
}
else {
// i == half-1
// Calculate the last "odd" prediction
int n_plus1 = new_n_plus1( vec[i-1], vec[i] );
predictVal = (int)((((float)vec[i] + (float)n_plus1)/2.0) + 0.5);
}
if (direction == forward) {
vec[j] = vec[j] - predictVal;
}
else if (direction == inverse) {
vec[j] = vec[j] + predictVal;
}
else {
printf("line::predict: bad direction value\n");
}
}
} // predict
/**
Update step of the linear interpolation wavelet
The predict phase works on the odd elements in the second
half of the array. The update phase works on the even
elements in the first half of the array. The update
phase attempts to preserve the average. After the update
phase is completed the average of the even elements should
be approximately the same as the average of the input data
set from the previous iteration. The result of the update
phase becomes the input for the next iteration.
In a Haar wavelet the average that replaces the even element is
calculated as the average of the even element and its neighboring
odd element (e.g., its odd neighbor before the split). In the
lifting scheme version of the Haar wavelet the odd element has
been overwritten by the difference between the odd element and
its even neighbor. In calculating the average (to replace the
even element) the value of the odd element can be recovered via
a simple algebraic manipulation.
In the line wavelet the odd element has been replaced by the
difference between the odd element and the mid-point of its two
even neighbors. Recovering the value of the odd element to
calculate the average is not as simple in this case.
The value that is added to the even element to preserve the
average is calculated by the equation shown below. This equation
is given in Wim Sweldens' journal articles and his tutorial
(Building Your Own Wavelets at Home) and in Ripples in
Mathematics. A somewhat more complete derivation of this
equation is provided in Ripples in Mathematics by A. Jensen
and A. la Cour-Harbo, Springer, 2001.
The equation used to calculate the average is shown below for a
given iteratin i. Note that the predict phase has already
completed, so the odd values belong to iteration i+1.
eveni+1,j = eveni,j op (oddi+1,k-1 + oddi+1,k)/4
This version of the line wavelet code implements an integer
version of linear interpolating wavelet. This versoin comes from
the paper Wavelet Transforms that Map Integers to Integers
by A.R. Calderbank, ingrid daubechies, wim weldens and Boon-Lock
Yeo, August 1996
This is the central reference that was used to develop this code.
Parts 1 and 2 of this paper are for the mathematicially
sophisticated (which is to say, they are not light reading).
However, for the implementer, part 3 and part 4 of this paper
provide excellent coverage of perfectly invertable wavelet
transforms that map integers to integers. In fact, part 3 of this
paper is worth reading in general for its discussion of the wavelet
Lifting Scheme.
The value added (or subtracted) from the eveni,j
(depending on whether the forward or inverse transform is being
calculated) is calculated from oddi+1,k-1 and
oddi+1,k from the predict step. This means that there
is missing value at the start of the set of odd elements (e.g., i
= 0, j == half). This missing value assumed to line on a line
with the first two odd elements.
Because interpolated values are used, the average is not
perfectly maintained.
*/
void update( T & vec, int N, transDirection direction )
{
int half = N >> 1;
for (int i = 0; i < half; i++) {
int j = i + half;
int val;
if (i == 0 && N == 2) {
val = (int)(((float)vec[j]/2.0) + 0.5);
}
else if (i == 0 && N > 2) {
int v_n_minus_1 = new_n_minus1( vec[j], vec[j+1] );
val = (int)((((float)v_n_minus_1 + (float)vec[j])/4.0) + 0.5);
}
else {
val = (int)((((float)vec[j-1] + (float)vec[j])/4.0) + 0.5);
}
if (direction == forward) {
vec[i] = vec[i] + val;
}
else if (direction == inverse) {
vec[i] = vec[i] - val;
}
else {
printf("update: bad direction value\n");
}
} // for
} // update
}; // line_int
#endif
packet/lossless/include/packcontainer_int.h 0100644 0000764 0000764 00000011651 07523321551 020131 0 ustar iank iank
#ifndef _PACKCONTAINER_INT_H_
#define _PACKCONTAINER_INT_H_
#include "packnode.h"
/**
The documentation in this file is formatted for doxygen
(see www.doxygen.org).
Copyright and Use
You may use this source code without limitation and without
fee as long as you include:
This software was written and is copyrighted by Ian Kaplan, Bear
Products International, www.bearcave.com, 2002.
This software is provided "as is", without any warranty or
claim as to its usefulness. Anyone who uses this source code
uses it at their own risk. Nor is any support provided by
Ian Kaplan and Bear Products International.
Please send any bug fixes or suggested source changes to:
iank@bearcave.com
@author Ian Kaplan
*/
/**
A container for use when calculating a packet wavelet tree.
By overriding the LHS and RHS versions of the [] operator,
packcontainer class allows a block of data to be treated as an array
while a wavelet transform step is being calculated. After the
wavelet transform step the data can be referenced as a left hand
side half (the result of the wavelet scaling function) and a right
hand side half (the result of the wavelet function). By allowing
the two halves of the array to be referenced, copying is avoided.
*/
class packcontainer_int {
private:
/** number of elements at this packet tree level */
size_t N;
/** left (low pass) half of the packcontainer_int data */
int* lhs;
/** right (high pass) half of the packnode data */
int* rhs;
private:
/** disallow the copy constructor */
packcontainer_int( const packcontainer_int &rhs ) {};
/** disallow default constructor */
packcontainer_int() {};
public:
/**
This version of the packcontainer object is used as a container
in calculating the forward wavelet packet transform for an integer
data set.
The packcontainer_int constructor is passed a pointer to a packnode
object. As the wavelet packet tree is being built this packnode
object contain the data from the previous level.
The length of the packnode object is N. The constructor will
dynamically allocate two vectors (lhs and rhs), each with N/2
elements. After the wavelet transform step has been calculated on
the packcontainer_int object, the lhs vector will contain the wavelet
scaling function, or low pass, result. The rhs vector will contain
the wavelet function, or high pass, result. These will be used
to construct two new packnode objects which will be children
of the object passed to the constructor.
*/
packcontainer_int( packnode* node )
{
assert( node != 0 );
N = node->length();
assert( N > 1 );
size_t half = N >> 1;
block_pool mem_pool;
size_t num_bytes = half * sizeof(int);
lhs = (int *)mem_pool.pool_alloc( num_bytes );
rhs = (int *)mem_pool.pool_alloc( num_bytes );
for (size_t i = 0; i < N; i++) {
(*this)[i] = (*node)[i];
}
} // packcontainer_int
/**
This version is used when calculating the inverse wavelet packet
transform. The constructor is passed the size of the container
(or at least the size that the container will be, when the left
and right hand size arrays are initialized). The lhs and rhs
arrays are then initialized and the inverse transform step
is calculated.
*/
packcontainer_int( size_t n )
{
N = n;
lhs = 0;
rhs = 0;
}
/**
This is a local implementation of new. When new
is applied to a packcontainer_int object, memory will
be allocated from a memory pool, rather than from
the system memory pool.
*/
void *operator new( unsigned int num_bytes )
{
block_pool mem_pool;
void *mem_addr = mem_pool.pool_alloc( num_bytes );
return mem_addr;
} // new
/** LHS [] operator */
int &operator[]( const size_t i )
{
assert( i < N );
size_t half = N >> 1;
if (i < half)
return lhs[i];
else {
return rhs[i-half];
}
}
/** RHS [] operator */
int operator[]( const size_t i ) const
{
assert( i < N );
size_t half = N >> 1;
if (i < half)
return lhs[i];
else {
return rhs[i-half];
}
}
/** return the left hand size array */
int* lhsData() { return lhs; }
/** return the right hand size array */
int* rhsData() { return rhs; }
/** set the left hand size array */
void lhsData(int* l) { lhs = l; }
/** set the right hand size array */
void rhsData(int* r) { rhs = r; }
/** return the length of the data in the packet
container. Note that this length is
the length of the rhs plus the length of
the lhs arrays.
*/
size_t length() { return N; }
}; // packcontainer_int
#endif
packet/lossless/include/packtree_base_int.h 0100644 0000764 0000764 00000003632 07523326251 020102 0 ustar iank iank
#ifndef _PACKTREE_BASE_INT_H_
#define _PACKTREE_BASE_INT_H_
/** \file
The documentation in this file is formatted for doxygen
(see www.doxygen.org).
Copyright and Use
You may use this source code without limitation and without
fee as long as you include:
This software was written and is copyrighted by Ian Kaplan, Bear
Products International, www.bearcave.com, 2002.
This software is provided "as is", without any warranty or
claim as to its usefulness. Anyone who uses this source code
uses it at their own risk. Nor is any support provided by
Ian Kaplan and Bear Products International.
Please send any bug fixes or suggested source changes to:
iank@bearcave.com
@author Ian Kaplan
*/
#include "packnode.h"
#include "liftbase.h"
#include "packcontainer_int.h"
/**
This is an integer version of the wavelet packet tree base
class. This class can be used for both compression and
time/frequency analysis of an integer time series. However,
this version was developed for lossless integer compression.
This class may be subclassed for classes that apply the
wavelet best basis algorithm or for time/frequency
wavelet packet trees.
*/
class packtree_base_int {
protected:
/** root of the wavelet packet tree */
packnode *root;
/** wavelet packet transform object */
liftbase *waveObj;
typedef enum { BadPrintKind,
printData,
printCost,
printBestBasis } printKind;
void breadthFirstPrint(printKind kind);
void newLevel( packnode* top, bool freqCalc, bool reverse );
public:
void pr();
/** get the root of the wavelet packet tree */
packnode *getRoot() { return root; }
}; // packtree_base_int
#endif
packet/lossless/include/packtree_int.h 0100644 0000764 0000764 00000005023 07523326600 017102 0 ustar iank iank
#ifndef _PACKTREE_INT_H_
#define _PACKTREE_INT_H_
#include "packtree_base_int.h"
#include "packdata_list.h"
#include "packcontainer.h"
#include "liftbase.h"
/** \file
The documentation in this file is formatted for doxygen
(see www.doxygen.org).
Copyright and Use
You may use this source code without limitation and without
fee as long as you include:
This software was written and is copyrighted by Ian Kaplan, Bear
Products International, www.bearcave.com, 2002.
This software is provided "as is", without any warranty or
claim as to its usefulness. Anyone who uses this source code
uses it at their own risk. Nor is any support provided by
Ian Kaplan and Bear Products International.
Please send any bug fixes or suggested source changes to:
iank@bearcave.com
@author Ian Kaplan
*/
/**
The packtree_int object constructs a wavelet packet tree
The constructor is passed a vector of integers, the length of the
vector (which must be a power of two) and a pointer to a wavelet
lifting scheme class that will be used in calculating the wavelet
transform step. If the vector passed to the constructor contains N
int vaues, the result of the constructor will be a wavelet packet
tree with log2(N) levels.
*/
class packtree_int : public packtree_base_int {
private:
/** Found original data marked as part of the best basis.
This means that the best basis function failed (or
that the original data is the most compact representation
relative to the cost function used). */
bool foundOriginalData;
/** found a best basis value in the wavelet packet tree */
bool foundBestBasisVal;
private:
/** disallow the copy constructor */
packtree_int( const packtree_int &rhs ) {};
/** disallow the default constructor */
packtree_int() {};
int bestBasisWalk( packnode *root );
void buildBestBasisList( packnode *root,
packdata_list &list );
void checkBestBasis( packnode *root );
void cleanTree(packnode *root, bool removeMark );
public:
packtree_int( const int *vec,
const size_t n,
liftbase *w );
/** destructor does nothing */
~packtree_int() {}
void prCost();
void prBestBasis();
void bestBasis();
bool bestBasisOK();
packdata_list getBestBasisList();
}; // packtree_int
#endif
packet/lossless/include/support.h 0100644 0000764 0000764 00000005573 07525261701 016162 0 ustar iank iank
#ifndef _SUPPORT_H_
#define _SUPPORT_H_
/** \file
The documentation in this file is formatted for doxygen
(see www.doxygen.org).
Copyright and Use
You may use this source code without limitation and without
fee as long as you include:
This software was written and is copyrighted by Ian Kaplan, Bear
Products International, www.bearcave.com, 2002.
This software is provided "as is", without any warranty or
claim as to its usefulness. Anyone who uses this source code
uses it at their own risk. Nor is any support provided by
Ian Kaplan and Bear Products International.
Please send any bug fixes or suggested source changes to:
iank@bearcave.com
@author Ian Kaplan
*/
/**
Various support functions for compression and bit width estimation.
These functions are the object oriented version of global functions.
They are all pure functions (e.g., no state in the function or
in the class) and static.
Bit width calculation
In compression applications an attempt is made to represent a data
set in a minimal number of bits. For these applications it must
also be possible to decompress the data. If the values in the
data set are stored in bit fields of different widths, information
about these fields must be stored as well.
Unlike compression algorithms, where decompresson must be taken
into account, this code simple calculates the sum of the bit
fields needed to represent a set of signed integer values.
This is useful for estimating how closely a wavelet transform
approximated a data set.
Double to integer conversion
The compression algorith uses integer to integer lossless
wavelet transforms. Since the financial data I am interested
in is in real form, it needs to be converted to integer.
This is done preserving three fractional digits.
The constructors and destructor are declared but not defined,
which will result in a link error if an instance of this class
is created (which is not the intent in the design).
*/
class support
{
private:
static size_t nearestPower2Width_( size_t val );
static int roundVal_( const double val );
public:
/** declare but do not define the constructor */
support();
/** declare but do not define the destructor */
~support();
/** declare but never define copy constructor */
support( const support &rhs );
static size_t valWidth( const int val );
size_t UnsignedValWidth( const size_t val );
static size_t vecWidth( const int *vec,
const size_t N );
static void roundToInt( int *intVec,
const double *realVec,
const size_t len );
static void decimalToInt( int *intVec,
const double *realVec,
const size_t len );
}; // support
#endif
packet/lossless/include/ts_trans_int.h 0100644 0000764 0000764 00000017755 07501537633 017166 0 ustar iank iank
#ifndef _TS_TRANS_INT_H_
#define _TS_TRANS_INT_H_
/** \file
The documentation in this file is formatted for doxygen
(see www.doxygen.org).
Copyright and Use
You may use this source code without limitation and without
fee as long as you include:
This software was written and is copyrighted by Ian Kaplan, Bear
Products International, www.bearcave.com, 2002.
This software is provided "as is", without any warranty or
claim as to its usefulness. Anyone who uses this source code
uses it at their own risk. Nor is any support provided by
Ian Kaplan and Bear Products International.
Please send any bug fixes or suggested source changes to:
iank@bearcave.com
@author Ian Kaplan
*/
#include "stdio.h"
#include "haar_int.h"
/**
An integer version of the TS transform (an extended S, or Haar
transform).
The TS transform is an extension of integer version of the Haar
trasform (which is sometimes referred the S-transform in the image
processing literature. The TS transform is an integer version of
the so called Cohen-Daubechies-Feaveau (3,1) transform. Here the
(3,1) refer to 3 vanishing moments for the wavelet function and 1
vanishing moment for the scaling function.
The equations for the lifting scheme version of the forward
TS transform are shown below. As with all lifting scheme
algorithms, the inverse transform exchanges addition and
subtraction operators.
The TS transform and the S transform are the same in the first two
steps (the first predict and update steps). An average
interpolation step is added to the S (Haar) transform.
d(1)1,i = s0,2i+1 - s0,2i
s1,i = s0,2i + floor( d(1)1,i/2 )
d1,i = d(1)1,i + floor((s1,i-1 - s1,i+1)/4.0 + 0.5)
This notation and the algorithm implemented here is taken directly
from Wavelet Transformst that Map Integers to Integers by
A.R. Calderbank, Ingrid Daugbechies, Wim Sweldens, and Boon-Lock
Yeo, 1996
The mathematical structure is reflected in the class structure.
Here the ts_trans_int class extends the haar_int class. The
haar_int class provide the predict() and update
functions. An additional predict2() function is added
that implements interpolation step.
Since an extra step has been added, the forwardStep and
inverseStep functions in the base class (liftbase)
are over-ridden by ts_trans_int.
A brief note on vanishing moments
This algorithm is commonly known as the CDF (3,1) wavelet
algorithm. As noted above, these numbers refer to the vanishing
moments. I have never found a definition of "vanishing moment" that
made intuitive sense to me, at least when applied to wavelet basis
functions.
As I understand the definition, a vanishing moment is a region over
which the integral is zero. So the area of sin(x) in the region
0..2Pi has an integral of zero, since the region between 0..Pi
results in a positive integral and the area between Pi..2Pi results
in the same integeral value, but with a negative sign. The sum of
these two regions is zero. If a wavelet were constructed from the
sine function, such that
x = {0..2Pi} : y = sin(x)
x = anything else : y = 0
This would be a wavelet with "compact support" (it is a defined for
a finite region, 0..2Pi) and one vanishing moment.
*/
// class ts_trans_int : public liftbase
class ts_trans_int : public haar_int
{
public:
/** the constructor does nothing */
ts_trans_int() {}
/** the destructor does nothing */
~ts_trans_int() {}
/** declare, but do not define the copy constructor */
ts_trans_int( const ts_trans_int &rhs );
private:
/**
Calculate the element s1,i+1.
The low pass half of the array is in the array index range
{0..half-1}. In the equation below, when i = half-1
a non-existent element at i+1 is needed. This function
calculates this element from s[half-2] and s[half-1].
d1,i = d(1)1,i + floor((s1,i-1 - s1,i+1)/4.0 + 0.5)
Here the non-existent element s[half] is assumed to lie on
the line from s[half-2] and s[half-1]. We pretend that
s[half-2] has the x-coordinate of 0 and s[half-1] has
an x-coordinate of 1. We then need to calculate the
y value at the x-coordinate 2.
The "two-point equation" for a line is used for this
calculation, where we are trying to find the value
of y, given
. y2 - y1
(y - y1) = -------- (x - x1)
. x2 - x1
Solving for y
. y2 - y1
y = -------- (x - x1) + y1
. x2 - x1
where
x1 = 0
x2 = 1
y1 = s[half-2]
y2 = s[half-1]
x = 2
Substituting in these values we get
y = 2*y2 - y1
*/
int new_n_plus1( int y1, int y2)
{
int y = 2 * y2 - y1;
return y;
}
/**
In the function new_n_plus1 a point beyond the end
of the low pass filter array is calculated. Here
a point beyond the beginning of the array is calculated.
x1 = 0
x2 = 1
y1 = s[0]
y2 = s[1]
x = -1
When these values are plugged into the
two point equation we get
y = 2 * y1 - y2
*/
int new_n_minus1( int y1, int y2)
{
int y = 2 * y1 - y2;
return y;
}
protected:
/**
The predict interpolation step.
Note that special cases exist at the start and end of the
array. A special case also exists for N = 2, where the
calculation is the same as the Haar wavelet (e.g., no
interpolation factor is added in).
*/
void predict2( int *& vec, int N, transDirection direction )
{
int half = N >> 1;
int predictVal;
for (int i = 0; i < half; i++) {
int j = i + half;
int y_n_plus1;
int y_n_minus1;
if (N == 2) {
y_n_minus1 = vec[0];
y_n_plus1 = vec[0];
}
else if (i == 0) {
y_n_minus1 = new_n_minus1( vec[0], vec[1] );
y_n_plus1 = vec[1];
}
else if (i < half-1) {
y_n_minus1 = vec[i-1];
y_n_plus1 = vec[i+1];
}
else { // i == half-1
y_n_minus1 = vec[i-1];
y_n_plus1 = new_n_plus1( vec[i-1], vec[i] );
}
predictVal = (int)( (((float)y_n_minus1 - (float)y_n_plus1)/4.0) + 0.5 );
if (direction == forward) {
vec[j] = vec[j] + predictVal;
}
else if (direction == inverse) {
vec[j] = vec[j] - predictVal;
}
else {
printf("haar_int::predict: bad direction value\n");
}
}
} // predict2
public:
/**
One TS transform forward step. This extends the S transform
(Haar integer transform) with the predict2 step.
*/
void forwardStep( int *& vec, const int n )
{
split( vec, n );
predict( vec, n, forward );
update( vec, n, forward );
predict2( vec, n, forward );
} // forwardStep
/**
One TS transform inverse step. This extends the S transform
(Haar integer transform) with the predict2 step.
*/
void inverseStep( int *& vec, const int n )
{
predict2( vec, n, inverse );
update( vec, n, inverse );
predict( vec, n, inverse );
merge( vec, n );
} // inverseStep
}; // ts_trans_int
#endif
packet/lossless/Makefile 0100644 0000764 0000764 00000005360 07525050265 014305 0 ustar iank iank
#
# Makefile for lossless wavelet compression software
#
# This makefile assumes that the Microsoft environment has
# been enhanced with the cygwin UNIX utilities for Win32.
# These are available from Red Hat (which purchased Cygnus).
#
# Command:
# nmake -f Makefile
#
DEBUG = -Zi
BROWSE = -FR
CFLAGS = $(BROWSE) $(DEBUG) -DWIN32 -Tp
OBJ = obj
EXE = .exe
TOPINC = ..\include
DOXYPATH = e:/doxygen/bin
INCLUDE = -Iinclude -I$(TOPINC) -I..\data\include -I "d:\Program Files\Microsoft Visual Studio\Vc98\Include"
PACKET_OBJ = blockpool.$(OBJ) packtree_base_int.$(OBJ) packtree_int.$(OBJ) invpacktree_int.$(OBJ) costwidth.$(OBJ) support.$(OBJ) yahooTS.$(OBJ)
all: compresstest
compresstest: compresstest$(EXE)
compresstest$(EXE): $(PACKET_OBJ) compresstest.$(OBJ)
$(CC) $(PACKET_OBJ) compresstest.$(OBJ) $(DEBUG) -o compresstest$(EXE)
#
# clean-up for Microsoft object and browser files and emacs temps
#
clean:
rm -f *.obj *.pdb *.sbr *.ilk *.exe
rm -f include/*~
rm -f src/*~
$(TOPINC)\costbase.h: $(TOPINC)\packnode.h
include\costwidth.h: $(TOPINC)\costbase.h
include\line_int.h: $(TOPINC)\liftbase.h
$(TOPINC)\packdata.h: $(TOPINC)\blockpool.h
include\invpacktree_int.h: $(TOPINC)\packdata.h
$(TOPINC)\packnode.h: $(TOPINC)\packdata.h
$(TOPINC)\packdata_list.h: $(TOPINC)\fifo_list.h $(TOPINC)\packdata.h
$(TOPINC)\packcontainer.h: $(TOPINC)\packnode.h
$(TOPINC)\queue.h: $(TOPINC)\fifo_list.h $(TOPINC)\packnode.h
include\packtree_base_int.h: $(TOPINC)\packnode.h $(TOPINC)\liftbase.h
include\packtree_int.h: include\packtree_base_int.h $(TOPINC)\packcontainer.h $(TOPINC)\liftbase.h
include\invpacktree_int.h: $(TOPINC)\liftbase.h $(TOPINC)\list.h $(TOPINC)\packcontainer.h $(TOPINC)\packdata.h $(TOPINC)\packdata_list.h
packtree_base_int.$(OBJ): src\$*.cpp include\$*.h
$(CC) -c $(INCLUDE) $(CFLAGS) src\$*.cpp
packtree_int.$(OBJ): src\$*.cpp include\$*.h $(TOPINC)\queue.h $(TOPINC)\packcontainer.h $(TOPINC)\blockpool.h
$(CC) -c $(INCLUDE) $(CFLAGS) src\$*.cpp
invpacktree_int.$(OBJ): src\$*.cpp include\$*.h $(TOPINC)\blockpool.h
$(CC) -c $(INCLUDE) $(CFLAGS) src\$*.cpp
blockpool.$(OBJ): ..\src\$*.cpp ..\include\$*.h
$(CC) -c $(INCLUDE) $(CFLAGS) ..\src\$*.cpp
costwidth.$(OBJ): src\$*.cpp include\$*.h
$(CC) -c $(INCLUDE) $(CFLAGS) src\$*.cpp
support.$(OBJ): src\$*.cpp include\$*.h
$(CC) -c $(INCLUDE) $(CFLAGS) src\$*.cpp
yahooTS.$(OBJ): ..\data\src\$*.cpp ..\data\include\$*.h
$(CC) -c $(INCLUDE) $(CFLAGS) ..\data\src\$*.cpp
compresstest.$(OBJ): src\$*.cpp include\line_int.h $(TOPINC)\packnode.h $(TOPINC)\packcontainer.h include\packtree_int.h include\invpacktree_int.h
$(CC) -c $(INCLUDE) $(CFLAGS) src\$*.cpp
packet/lossless/src/ 0040755 0000764 0000764 00000000000 07525263472 013441 5 ustar iank iank packet/lossless/src/compresstest.cpp 0100644 0000764 0000764 00000017241 07525263247 016702 0 ustar iank iank
/** \file
This file contains code to test various lossless compression
algorithms on financial time series (e.g., stock close price).
The documentation in this file is formatted for doxygen
(see www.doxyeng.org).
Copyright and Use
You may use this source code without limitation and without
fee as long as you include:
This software was written and is copyrighted by Ian Kaplan, Bear
Products International, www.bearcave.com, 2002.
This software is provided "as is", without any warranty or
claim as to its usefulness. Anyone who uses this source code
uses it at their own risk. Nor is any support provided by
Ian Kaplan and Bear Products International.
Please send any bug fixes or suggested source changes to:
iank@bearcave.com
@author Ian Kaplan
*/
#include
#include "invpacktree_int.h"
#include "packtree_int.h"
#include "support.h"
#include "delta.h"
#include "haar_int.h"
#include "ts_trans_int.h"
#include "line_int.h"
#include "costwidth.h"
#include "yahooTS.h"
/**
Make a copy of an integer array. Note that the function
allocates memory, but does not deallocate it.
*/
int *copy( int *intVec, const size_t N )
{
int *newVec = new int[ N ];
for (size_t i = 0; i < N; i++)
{
newVec[i] = intVec[i];
}
return newVec;
} // copy
/**
Compare two integer arrays of size N. Return true if they are
equal, false otherwise.
*/
bool compare( const int *v1, const int *v2, const size_t N )
{
bool rslt = true;
for (size_t i = 0; i < N; i++)
{
if (v1[i] != v2[i])
{
rslt = false;
break;
}
}
return rslt;
} // compare
/**
Calculate the number of bits needed to represent the result
of the wavelet packet transform.
The result of this version of the wavelet packet transform
is a list of vectors, where each vector represents a best basis
fit for a particular region of the signal.
The result returned here is an obvious lower bound, since in
practice a vector would consist of a fixed set of values. Also
the length of the vector would have to be included, along
with the width the the length values.
*/
size_t calcPacketWidth( packdata_list &bestBasisList )
{
packdata_list::handle h;
size_t totalWidth = 0;
for (h = bestBasisList.first(); h != 0; h = bestBasisList.next( h )) {
packdata *node = bestBasisList.get_item( h );
const size_t len = node->length();
const int *vec = node->getData();
int nodeWidth = support::vecWidth( vec, len );
totalWidth += nodeWidth;
}
return totalWidth;
} // calcPacketWidth
/**
Wavelet packet calculation
The function returns the minimal total number of bits needed
to represent the data when compressed using the wavelet
packet algorithm.
*/
size_t packet_calc( const int *intVec,
const int *copyVec,
const int N )
{
// The "line" wavelet transform (e.g., line with slope)
line_int line;
// calculate the wavelet packet tree, using the line wavelet transform
packtree_int tree( intVec, N, &line );
// get the root of the wavelet packet transform tree
packnode *treeRoot = tree.getRoot();
// Assign a cost on the basis of bit width
costwidth cost( treeRoot );
// Calculate the "best basis" function from the tree
tree.bestBasis();
// Check that the best basis function succeeded. That is,
// that the best basis function does not include the
// original data.
if (! tree.bestBasisOK()) {
printf("Best basis calculation failed\n");
}
// Get the best basis list. This will be a list of
// nodes consisting of the best basis set. This set is
// obtained by traversing the tree, top down, left to
// right.
packdata_list bestBasis = tree.getBestBasisList();
// Sum the cost values (width for each node in the best basis
// list
size_t width = calcPacketWidth( bestBasis );
invpacktree_int invtree( bestBasis, &line );
const int *invRslt = invtree.getData();
bool isEqual = compare( invRslt, copyVec, N );
if (! isEqual)
printf("Wavelet packet inverse is wrong\n");
return width;
} // packet_calc
/**
Calculate the number of bits needed to represent the data after
"delta" compression. Delta compression stores the difference
between value si-1 and si.
*/
size_t delta_calc( int *intVec,
const int *copyVec,
const int N )
{
delta_trans delta;
delta.forward( intVec, N );
const size_t deltaWidth = support::vecWidth( intVec, N );
delta.inverse( intVec, N);
bool isEqual = compare( intVec, copyVec, N);
if (! isEqual)
printf("Delta compression inverse failed\n");
return deltaWidth;
} // delta_calc
/**
Calculate the number of bits needed to represent the data
after wavelet compression wavelet compression.
\param intVec A pointer to an array of integers
\param copyVec A pointer to a copy of the integer array.
This is used to verify that the inverse(forward(intVec))
transform yields the original data.
\param N the size of intVec and copyVec
\param w a pointer to a wavelet class
*/
size_t wave_calc( int *intVec,
const int *copyVec,
const int N,
liftbase *w )
{
w->forwardTrans( intVec, N );
const size_t waveWidth = support::vecWidth( intVec, N );
w->inverseTrans( intVec, N );
bool isEqual = compare( intVec, copyVec, N);
if (! isEqual)
printf("Line wavelet inverse is wrong\n");
return waveWidth;
} // wave_calc
/**
Read in a set of equity time series file. Calculate the
number of bits needed to represent the data without
compression and after wavelet and wavelet packet
compression.
*/
int main()
{
const char *files[] = { "aa", // Alcoa Aluminium
"amat", // Applied Materials
"ba", // Boeing
"cof", // Capital One
"ge", // General Electric
"ibm", // IBM Corp.
"intc", // Intel
"mmm", // 3M
"mrk", // Merck
"wmt", // Wal-Mart
0 // The null pointer
};
const size_t N = 512;
double realVec[ N ];
int intVec[ N ];
// an instance of yahooTS with the path to the data directory
const char *dataDirPath = "..\\data\\equities\\";
yahooTS ts( dataDirPath );
printf("Equity Uncompressed delta Haar line TS wavelet packet (line)\n");
for (size_t i = 0; files[i] != 0; i++) {
size_t n = N;
if (! ts.getTS( files[i], realVec, n, yahooTS::Close )) {
break;
}
if (n != N) {
printf("Error: %d out of %d data elements read\n", n, N );
break;
}
support::decimalToInt( intVec, realVec, N );
int *copyVec = copy( intVec, N );
const size_t beforeWidth = support::vecWidth( intVec, N );
const size_t deltaWidth = delta_calc( intVec, copyVec, N );
haar_int haar;
const size_t haarWidth = wave_calc( intVec, copyVec, N, &haar );
line_int line;
const size_t lineWidth = wave_calc( intVec, copyVec, N, &line );
ts_trans_int ts_trans;
const size_t tsTransWidth = wave_calc( intVec,
copyVec,
N,
&ts_trans);
const size_t packetWidth = packet_calc( intVec, copyVec, N );
printf(" %4s %4d %4d %4d %4d %d %4d\n",
files[i], beforeWidth, deltaWidth, haarWidth, lineWidth, tsTransWidth, packetWidth );
}
return 0;
}
packet/lossless/src/costwidth.cpp 0100644 0000764 0000764 00000002023 07523330015 016131 0 ustar iank iank
/** \file
The documentation in this file is formatted for doxygen
(see www.doxygen.org).
Copyright and Use
You may use this source code without limitation and without
fee as long as you include:
This software was written and is copyrighted by Ian Kaplan, Bear
Products International, www.bearcave.com, 2002.
This software is provided "as is", without any warranty or
claim as to its usefulness. Anyone who uses this source code
uses it at their own risk. Nor is any support provided by
Ian Kaplan and Bear Products International.
Please send any bug fixes or suggested source changes to:
iank@bearcave.com
@author Ian Kaplan
*/
#include "costwidth.h"
#include "support.h"
int
costwidth::costCalc( packnode *root )
{
assert( root != 0 );
size_t N = root->length();
const int *a = root->getData();
int width = support::vecWidth( a, N );
return width;
} // costCalc
packet/lossless/src/invpacktree_int.cpp 0100644 0000764 0000764 00000015403 07523334337 017327 0 ustar iank iank
#include
#include
#include "blockpool.h"
#include "invpacktree_int.h"
/** \file
For additional information on this implementation of the inverse
wavelet packet transform see
http://www.bearcave.com/misl/misl_tech/wavelets/packet/index.html
The documentation in this file is formatted for doxygen
(see www.doxyeng.org).
Copyright and Use
You may use this source code without limitation and without
fee as long as you include:
This software was written and is copyrighted by Ian Kaplan, Bear
Products International, www.bearcave.com, 2002.
This software is provided "as is", without any warranty or
claim as to its usefulness. Anyone who uses this source code
uses it at their own risk. Nor is any support provided by
Ian Kaplan and Bear Products International.
Please send any bug fixes or suggested source changes to:
iank@bearcave.com
@author Ian Kaplan
*/
/**
Allocate a new level in the wavelet packet tree that is being
constructed. The new "level" is built from a packcontainer_int object.
The left child will be the elem argument. The length of the data in
the new level will be twice that of lower level (of which elem is
half).
*/
void invpacktree_int::new_level( packdata *elem )
{
size_t half = elem->length();
size_t n = half * 2;
packcontainer_int *container = new packcontainer_int( n );
container->lhsData( (int *)elem->getData() );
stack.add( container );
} // new_level
/**
At this point the Top Of Stack (TOS) packcontainer_int object should
have both a right and a left array. Calculate an inverse wavelet
transform step on the packcontainer_int object. The packcontainer_int
object allows the right and left hand side arrays to be treated as
one array
If the current top of stack is twice the size of the inverse wavelet
transform step result, the result becomes the right hand size of
top of stack packcontainer_int and reduce is called recursively.
If the TOS is empty or it is not twice the size of the inverse
transform result, a new packcontainer_int will be pushed on the stack.
The left hand size will be the transform result.
*/
void invpacktree_int::reduce()
{
LIST::handle h;
h = stack.first();
packcontainer_int *tos = stack.get_item( h );
assert( tos->lhsData() != 0 && tos->rhsData() != 0 );
/** Remove the linked list element that used to contain tos
(e.g., pop the linked list element off). Note that this
leaves the tos object unchanged (e.g., it does not delete
it).
*/
stack.remove();
size_t n = tos->length();
// calculate the inverse wavelet transform step
waveObj->inverseStep( (*tos), n );
// copy the result of the inverse wavelet transform
// into a new data array.
block_pool mem_pool;
int *vec = (int *)mem_pool.pool_alloc( n * sizeof( int ) );
for (int i = 0; i < n; i++) {
vec[i] = (*tos)[i];
}
if (stack.first() != 0) {
h = stack.first();
packcontainer_int *tos = stack.get_item( h );
if (tos->length() == n*2) {
tos->rhsData( vec );
reduce();
}
else {
assert( tos->length() > n*2 );
packcontainer_int *container = new packcontainer_int( n*2 );
container->lhsData( vec );
stack.add( container );
} // else
}
else {
// the stack is empty
packcontainer_int *container = new packcontainer_int( n*2 );
container->lhsData( vec );
stack.add( container );
}
} // reduce
/**
Add a packdata element to the inverse wavelet packet transform
calculation.
A packdata element contains the wavelet packet data items from
one level of the wavelet packet transform tree.
The inverse wavelet packet transform calculation uses a stack. Each
level in the stack consists of a packcontainer_int object. This object
consists of left and right hand size arrays. When both these arrays
are present, they can be treated as one contiguous array (courtesy
of C++ operator overloading). The size of the packcontainer_int object
is twice the size of its left and right hand size arrays.
If the stack is empty, a new level is created. A packetcontainer
object can be viewed as a binary tree element. The left hand
size is filled in first.
If the stack is not empty and the packcontainer_int object on the
top of stack (TOS) is twice the size of the elem argument,
then the array contained in the elem argument is added to
the TOS element and reduce is called to calculate a step
of the inverse wavelet transform.
If the TOS element is greater than twice the size of elem
then a new level is added.
*/
void invpacktree_int::add_elem( packdata *elem )
{
assert( elem != 0 );
if (stack.first() == 0) {
new_level( elem );
}
else {
size_t half = elem->length();
size_t n = half * 2;
LIST::handle h;
h = stack.first();
packcontainer_int *tos = stack.get_item( h );
if (tos->length() == n) {
assert( tos->rhsData() == 0);
tos->rhsData( (int *)elem->getData() );
reduce();
}
else if (tos->length() > n) {
new_level( elem );
}
else {
printf("add_elem: the size of the TOS elem is wrong\n");
}
} // else
} // add_elem
/**
This constructor calculates the inverse wavelet packet transform.
The constructor is passed a packdata_list object, which contains
the "best basis" result of the wavelet packet transform. The
liftbase template argument is a pointer to a wavelet transform
function. This wavelet transform must be the same function that
was used to calculate the packet transform.
*/
invpacktree_int::invpacktree_int( packdata_list &list,
liftbase *w )
{
data = 0;
N = 0;
waveObj = w;
// Traverse the "best basis" list and calculate the inverse
// wavelet packet transform.
packdata_list::handle h;
for (h = list.first(); h != 0; h = list.next( h )) {
packdata *elem = list.get_item( h );
add_elem( elem );
} // for
LIST::handle tosHandle;
tosHandle = stack.first();
packcontainer_int *tos = stack.get_item( tosHandle );
if (tos != 0) {
size_t len = tos->length();
N = len/2;
data = tos->lhsData();
stack.remove();
}
} // invpacktree_int
/**
print the result of the inverse wavelet packet transform
*/
void invpacktree_int::pr()
{
if (data != 0) {
for (int i = 0; i < N; i++) {
printf("%7.4f ", data[i] );
}
printf("\n");
}
} // pr
packet/lossless/src/packtree_base_int.cpp 0100644 0000764 0000764 00000013140 07523321622 017571 0 ustar iank iank
#include "queue.h"
#include "packnode.h"
#include "packcontainer_int.h"
#include "packtree_base_int.h"
/**
Add a new level to the wavelet packet tree. Wavelet packet trees
are data structures that support a variety of applications.
If the reverse argument is true, the locations of the high pass and
low pass filter results in the wavelet calculation will be
reversed. This is used in building a wavelet packet tree for
frequency analysis.
The top packnode object contains data from the previous
level. If this is the first level in the tree, top will
contain the input data set.
A packcontainer_int object is created with the top data. The
packcontainer_int constructor will allocate new storage and copy the
data into this storage. If the packcontainer_int object consists of N
elements, then there will be N/2 elements on the left and N/2 or the
right. The object allows the lhs and rhs data to be accessed as one
array.
A wavelet transform step is calculated on the N element data set.
This results in N/2 values from the wavelet scaling function (low
pass function) and N/2 values from the wavelet function (high pass
function). These values are used to create two new packnode objects
which become children of top. Sub-trees for the new packnode
objects are recursively calculated.
As the tree is constructed, the leaves of the tree are marked with a
boolean flag in preparation for calculating the "best basis"
representation for the data. See the algorithm outlined in section
8.2.2 of "Ripples in Mathematics" by Jensen and la Cour-Harbo
The wavelet packet tree form that is used for wavelet packet
frequency analysis is described in section 9.3 (figure 9.14)
ofg "Ripples in Mathematics".
*/
void packtree_base_int::newLevel( packnode* top,
bool freqCalc,
bool reverse )
{
if (top != 0) {
const size_t len = top->length();
if (len > 1) {
// Create a new wavelet packet container for use in
// calculating the wavelet transform. Note that the
// container is only used locally.
packcontainer_int container( top );
if (reverse) {
// Calculate the reverse foward wavelet transform step,
// where the high pass result is stored in the upper half
// of the container and the low pass result is stored
// in the lower half of the container.
waveObj->forwardStepRev( container, len );
}
else {
// Calculate the foward wavelet transform step, where
// the high pass result is stored in the upper half
// of the container and the low pass result is stored
// in the lower half of the container.
waveObj->forwardStep( container, len );
}
packnode *lhs = new packnode(container.lhsData(),
len/2,
packnode::LowPass );
packnode *rhs = new packnode(container.rhsData(),
len/2,
packnode::HighPass );
// set the "mark" in the top node to false and
// mark the two children to true.
top->mark( false );
lhs->mark( true );
rhs->mark( true );
top->lhsChild( lhs );
top->rhsChild( rhs );
// The transform on the left hand side always uses
// the standard order (e.g., low pass filter result
// goes in the lower half, high pass goes in the
// upper half of the container).
newLevel( lhs, freqCalc, false );
if (freqCalc) {
// wavelet packet frequency analysis reverses the
// storage locations for the filter results in the
// right hand child
newLevel( rhs, freqCalc, true );
}
else { // freq == false
// use standard filter location
newLevel( rhs, freqCalc, false );
}
}
}
} // newLevel
/**
Print the wavelet packet tree, breadth first (this is also
sometimes called a level traversal).
The breadth first traversal uses a queue. The root of the tree is
initially inserted into the queue. The function operates by
deleting the node at teh front of the queue, printing the data
associated with that node and adding the node's left and right
children to the queue. Since a node's children are at the next
lower level, and we add the left child before the right child,
the function prints a tree level from left to right.
The above paragraph paraphrases Chapter 5, Level Order Traversal
of Fundamentals of Data Structures in C by Horowitz, Sahni
and Anderson-Freed, 1993.
*/
void packtree_base_int::breadthFirstPrint(const printKind kind)
{
queue Q;
Q.addQueue( root, 0 );
while (! Q.queueEmpty() ) {
packnode *node = Q.queueStart()->node;
size_t indent = Q.queueStart()->indent;
Q.deleteStart();
if (indent > 0) {
// print 'indent' spaces
printf("%*c", indent, ' ');
}
switch (kind) {
case printData: node->pr();
break;
case printCost: node->prCost();
break;
case printBestBasis: node->prBestBasis();
break;
default:
assert( false );
break;
} // switch
packnode *lhs = node->lhsChild();
packnode *rhs = node->rhsChild();
if (lhs != 0) {
Q.addQueue( lhs, indent + 2 );
}
if (rhs != 0) {
Q.addQueue( rhs, indent + 2 );
}
}
} // packtree_base_int::breadthFirstPrint
/**
Print the wavelet packet tree data and wavelet transform
result to standard out.
*/
void packtree_base_int::pr()
{
if (root != 0) {
breadthFirstPrint(printData);
}
} // pr
packet/lossless/src/packtree_int.cpp 0100644 0000764 0000764 00000017152 07523326656 016622 0 ustar iank iank
/** \file
The documentation in this file is formatted for doxygen
(see www.doxyeng.org).
Copyright and Use
You may use this source code without limitation and without
fee as long as you include:
This software was written and is copyrighted by Ian Kaplan, Bear
Products International, www.bearcave.com, 2002.
This software is provided "as is", without any warranty or
claim as to its usefulness. Anyone who uses this source code
uses it at their own risk. Nor is any support provided by
Ian Kaplan and Bear Products International.
Please send any bug fixes or suggested source changes to:
iank@bearcave.com
@author Ian Kaplan
*/
#include
#include
#include "packcontainer_int.h"
#include "packtree_int.h"
/**
Construct a wavelet packet tree from a vector of integer
values. The size of the vector, which must be a power
of two, is passed in N. A wavelet Lifting Scheme object
is passed in w. This object is used to calculate
the wavelet transform step which is applied at each level
(where level > 0) of the wavelet packet tree.
The first level (level 0) of the wavelet packet tree contains
the original data set.
\arg vec An array of integer values on which the wavelet packet
transform is calculated.
\arg N The number of elements in the input array
\arg w A pointer to the the wavelet transform object to use
in calculating the wavelet packet transform.
*/
packtree_int::packtree_int( const int *vec,
const size_t N,
liftbase *w )
{
waveObj = w;
block_pool mem_pool;
int *vecCopy = (int *)mem_pool.pool_alloc( N * sizeof( int ) );
for (int i = 0; i < N; i++) {
vecCopy[i] = vec[i];
}
root = new packnode( vecCopy, N, packnode::OriginalData );
root->mark( true );
newLevel( root, false, false );
} // packtree_int
/**
The best basis algorithm selects the nodes nearest the tree root for
the best basis set. These nodes are "marked" true with a boolean
flag.
The best basis algorithm outlined in Ripples in Mathematics
sets any marks in nodes that are below a marked node in the tree
(nodes which are in a subtree of a marked node) to false. However,
this is unnecessary when the best basis set is collected, since the
algorithm uses top down tree traversal and stops at marked node.
A problem does arise when the entire tree is printed to show the
nodes that are marked as part of the best basis set. In this case
the result may appear wrong, since child nodes of a best basis node
are marked. This function does a top down traversal setting the
"marks" on these child nodes to false.
*/
void packtree_int::cleanTree(packnode *top, bool removeMark )
{
if (top != 0) {
if (removeMark) {
if (top->mark()) {
top->mark( false );
}
}
else {
// if a mark is found, then set all the "marks" below that
// point to false (e.g., remove the marks).
if (top->mark()) {
removeMark = true;
}
}
cleanTree( top->lhsChild(), removeMark );
cleanTree( top->rhsChild(), removeMark );
}
} // cleanTree
/**
Print the wavelet packet tree cost values in breadth first
order.
*/
void packtree_int::prCost()
{
if (root != 0) {
breadthFirstPrint(printCost);
}
}
/**
Print the wavelet packet tree, showing the nodes
that have been selected by the "best basis" algorithm.
*/
void packtree_int::prBestBasis()
{
if (root != 0) {
cleanTree( root, false );
breadthFirstPrint(printBestBasis);
}
} // prBestBasis
/**
Walk the wavelet packet tree and apply the "best basis" algorithm
described in Chapter 8 of Ripples in Mathematics. Nodes that
are "marked" become part of the best basis, which is a minimal
representation of the data in terms of the cost function.
*/
int packtree_int::bestBasisWalk( packnode *top )
{
int cost = 0;
if (top != 0) {
packnode *lhs = top->lhsChild();
packnode *rhs = top->rhsChild();
if (lhs == 0 && rhs == 0) { // we've reached a leaf
cost = top->cost();
}
else if (lhs != 0 && rhs != 0) {
int lhsCost = bestBasisWalk( lhs );
int rhsCost = bestBasisWalk( rhs );
int v1 = top->cost();
int v2 = lhsCost + rhsCost;
if (v1 <= v2) {
top->mark( true );
lhs->mark( false );
rhs->mark( false );
}
else { // v1 > v2
top->cost( v2 );
}
cost = top->cost();
}
else {
// The tree does not seem to be a full binary tree
// Something has gone badly wrong.
assert( false );
}
}
return cost;
} // bestBasicWalk
/**
Calculate the wavelet packet "best basis"
*/
void packtree_int::bestBasis()
{
bestBasisWalk( root );
} // bestBasis
/**
Recursively traverse the wavelet packet tree and check that the best
basis result is correct. That is, that the best basis has been
calculated and that it does not include the orignal data set.
The best basis includes the original data set with the packnode
ofr the original data set is marked.
This algorithm makes use of two class global variables. There may
be a purely recursive, way to do this without these global class
variables, but these variables make the algorithm much easier.
The variables are initialized by the calling function
bestBasisOK().
*/
void packtree_int::checkBestBasis( packnode *top )
{
if (top != 0) {
if (top->mark()) {
foundBestBasisVal = true;
if (top->getKind() == packdata::OriginalData) {
foundOriginalData = true;
}
}
if (!foundOriginalData) {
checkBestBasis( top->lhsChild() );
}
if (!foundOriginalData) {
checkBestBasis( top->rhsChild() );
}
}
} // checkBestBasis
/**
Return true is be best basis has been calculated properly,
return false if the best basis has not been calculated or
it was not calculated properly.
The best basis is calculated in reference to a particular
cost function. A particular cost function will not always
result in a data set which differs from the original data
set. If this is the case, the "best basis" result will
include the original data.
*/
bool packtree_int::bestBasisOK()
{
foundOriginalData = false;
foundBestBasisVal = false;
checkBestBasis( root );
bool rslt = (foundBestBasisVal && (!foundOriginalData));
return rslt;
} // bestBasicOK
/**
Traverse the tree from the top down and add the best basis
nodes to the best basis list.
Note that the list object is simply a package for a scalar
value, the pointer to the head of the list. So it can be
passed by value without incurring a cost greater than
passing a pointer (e.g., pass by reference).
*/
void packtree_int::buildBestBasisList( packnode *top,
packdata_list &list )
{
if (top != 0) {
if (top->mark()) {
list.add( top );
}
else {
buildBestBasisList( top->lhsChild(), list );
buildBestBasisList( top->rhsChild(), list );
}
}
} // buildBestBasisList
/**
Return a list consisting of the best basis packdata values.
*/
packdata_list packtree_int::getBestBasisList()
{
packdata_list list;
buildBestBasisList( root, list );
return list;
} // getBestBasisList
packet/lossless/src/support.cpp 0100644 0000764 0000764 00000011152 07525263210 015644 0 ustar iank iank
#include "support.h"
/**
Represent a real value with three fractional digits in integer
form. To accomplish this the real number is rounded and then
multiplied by 1000.0.
The rounding used is so proper rounding:
-
If the digit in the fourth
place is greater than 5, the digit in the third place is incremented
by one.
-
If the digit in the fourth place is less than 5, the third digit
is left unchanged.
-
If the digit in the fourth place is 5 and the digit in the third
place is odd, the third digit is incremented.
Proper rounding is unbiased (that is, in a random sample of numbers
the same number of values are rounded up as are rounded down).
Some examples:
12.4567 is rounded to 12.457 and converted to the integer
12457.
42.1234 is rounded to 42.123 and converted to the integer
42123.
127.1235 is rounded to 127.124 and converted to the integer
127124.
While this function does what I intended, it does seem to do it
in a lot of operations. There is probably a faster way to do this
by directly manipulating IEEE floating point. While this might
be faster, it is more complex.
*/
int support::roundVal_( const double val )
{
int intPart = static_cast(val);
double fracPart = val - intPart;
int threeDigits = fracPart * 1000.0;
int fourDigits = fracPart * 10000.0;
int forthDigit = fourDigits % 10;
int thirdDigit = threeDigits % 10;
double roundVal = 0.001;
if (forthDigit < 5) {
roundVal = 0.0;
}
else if (forthDigit == 5) {
if ((thirdDigit & 0x1) == 0) {
roundVal = 0.0;
}
}
double newVal = val + roundVal;
double intRslt = newVal * 1000.0;
return intRslt;
} // roundVal_
/**
Round an array of doubles to three decimal places and multiply
by 1000, resulting in an integer that reflects the double
value.
For equity time series this kind of rounding is useful for
"old style" pre-decimalization time series where the fractional
values represent increments by 1/16. This kind of rounding is
also useful for adjusted time series (e.g., the time series
is adjusted by a split and/or dividends are reinvested).
Since this function results in an integer vector with integer
values reflecting three decimal places, this function should
not be used for decimalized values (which only have two
decimal places).
*/
void support::roundToInt( int *intVec,
const double *realVec,
const size_t len )
{
if (intVec != 0 && realVec != 0) {
for (size_t i = 0; i < len; i++) {
intVec[i] = roundVal_( realVec[i] );
}
}
} // roundToInt
/**
Convert a decimalized array of doubles (e.g., where there
are two base ten fractional digits (e.g., 6.02, 3.14, 1.15)
into integer form. This is done my multiplying by 100 and
truncating. No rounding is done, since this function
assumes that here are only two significant fractional digits.
The roundToInt function can be used for rounding.
*/
void support::decimalToInt( int *intVec,
const double *realVec,
const size_t len )
{
if (intVec != 0 && realVec != 0) {
for (size_t i = 0; i < len; i++) {
intVec[i] = (int)(realVec[i] * 100.0);
}
}
} // roundToInt
size_t support::nearestPower2Width_( size_t val )
{
size_t width = 0;
if (val > 0) {
width = 1;
size_t power = 1;
while (power < val && width < 32) {
power = power << 1;
width++;
}
}
return width;
} // nearestPower2Width_
/**
Calculate the number of bits needed to represent an integer
value. A sign bit is added, so that positive numbers always
have a leading zero and negative numbers have a leading one.
*/
size_t support::valWidth( const int val )
{
size_t wholeNum = (val < 0) ? -val : val;
size_t width = 1 + nearestPower2Width_( wholeNum );
return width;
} // valWidth
/**
Calculate the number of bits needed to represent an unsigned
value.
*/
size_t support::UnsignedValWidth( const size_t val )
{
size_t width = nearestPower2Width_( val );
return width;
} // valWidth
/**
Calculate the minimum number of bits needed to represent the
values in an integer vector.
*/
size_t support::vecWidth( const int *vec, const size_t N )
{
size_t totalWidth = 0;
if (vec != 0) {
for (size_t i = 0; i < N; i++) {
totalWidth += valWidth( vec[i] );
}
}
return totalWidth;
} // vecWidth( int *)
packet/MainPage 0100644 0000764 0000764 00000027052 07525057655 012416 0 ustar iank iank
/*! \mainpage
Introduction
This the main page for the doxygen
documentation for the following C++ source code that implements:
-
Wavelet packet algorithm
-
Inverse wavelet packet algorithm
-
Time/Frequency analysis version of the wavelet packet algorithm
-
Floating point (e.g., double) versoins of the Daubechies D4, Haar
and "line" wavelet algorithms.
-
Lossless integer wavelet compression algorithms (Haar, TS Transform
and "line" wavelets).
-
Integer wavelet packet transform (for lossless compression)
-
Support code to read in Yahoo historical equity data files.
The C++ wavelet algorithms published here are related to a
set of Java algorithms also published on bearcave.com (see
Wavelets
in Java).
C++ was used to implement the wavelet packet algorithms because
C++ provides features like operator overloading and templates
which proved useful in creating reusable wavelet packet code.
The associated Web pages can be found on Bearcave.com at:
The web pages above provide an overview of the wavelet packet
transform. These web pages are meant to be read along with this
source code, which is extensively documented.
The wavelet packet transform has a variety of applications, but one
of the most common involves data compression. A cost function is
associated with the wavelet packet transform which will find the
best wavelet "basis" for various regions of the data set. The "best
basis" is the wavelet scale that best fits that region. The
determination of "best fit" is made relative to a given wavelet and
cost function pair.
Assuming that the wavelet function can produce an approxiation of
the data set, the end result of both the floating point and integer
versions of the wavelet packet transform will be the data set mapped
into a set of smaller values. In the case of lossy compression
small values can be set to zero, allowing coding compression (e.g.,
Huffman or "run length" encoding).
The integer version of the wavelet packet algorithm (which uses
lossless integer to integer wavelet transforms) is useful for
lossless compression. In this case, compression is achieved by the
fact that the data set can be mapped to smaller values which
can be represented in fewer bits.
The source code included here also includes a modified
wavelet packet algorithm, which is used for time/requency
analysis.
Given the complexity of wavelets, the web pages above (and the
commented source code) fall short in completeness of explanation.
The material I've written should be read along with other sources
(e.g., books and articles on wavelets). No single reference will
cover all the facets of wavelet mathematics and application. There
is a vast literature on wavelets. Many articles and books are
written by mathematicians, either for other mathematicians or for
students of mathematics. One of the few books I've found that leans
less toward theory and more toward application is Ripples in
Mathematics by Jensen and la-Cour Harbo, Springer Verlag, 2001.
The wavelet packet algorithms implemented here are based on Chapter
8 and 9 of this book.
Download
You can download the source code, in tar format, for the wavelet
packet transform
here. The source include Makefiles for both Windows NT (nmake) and
UNIX.
This source code contains both the wavelet packet transform code,
the modified wavelet packet transform for wavelet time
frequency analysis, lossless wavelet compression code and
support to read Yahoo equity time series files.
Trees vs. Tables
I have chosen to represent the result of the wavelet packet
transform as a dynamically allocated tree. The same data
can be represented as a table, where each table row corresponds
to a horizontal slice through the tree (e.g., the level basis).
I think that the tree representation is more flexible and natural.
Dynamic data structures are easily implemented in languages like C++
and Java. This may not be true of MATLAB code.
Overview of the Source Code
The source code published here can be divided into three catagories:
-
Floating point wavelet packet and modified wavelet packet code.
The wavelet packet code uses floating point wavelet algorithms
and floating point cost functions (e.g., Shannon and threshold
cost functions).
-
Lossless wavelet packet compression code. This code is supported
by integer to integer wavelet transforms.
-
Memory management. The wavelet packet algorithm uses dynamic
allocation (e.g., new). To make deallocation fast,
memory is allocated from a memory pool. This allows all
allocated memory to be deallocated at once. Pool allocation
also simplifies the code, since a single deallocation call
is made when the data structures are no longer needed.
-
Data input support for the lossless wavelet compression
algorithm. My motivation for developing the lossless wavelet
compression code is the application of this algorithm to financial
time series (e.g., stock close prices). In this case these are
historical data sets downloaded from Yahoo.
Data structure templates
The wavelet packet algorithm makes use of a several
general data structure templates:
-
FIFO_LIST template. This template supports a list with
both a list head and a tail tail pointer. Items are added to the
end of the list. When read from the front, the items will
be read in first-in, first-out order. This is used to build a
"queue" which is used in calculating the inverse wavelet packet
transform and in constructing the wavelet packet data list that
results from the wavelet packet transform.
-
LIST template. This supports a simple linked list
with a single list head pointer.
-
GrowableArray template. This is a generic growable array
that is by the modified wavelet packet algorithm for time
frequency analysis.
Wavelet Algorithm Templates
The wavelet algorithms used by the wavelet packet algorithms are based
on the wavelet lifting scheme, where thare are one or more predict and
update stages. The liftbase template provides a common base
class and can be instantiated for the array type and the array element
type. The wavelet algorithms themselfs are also templatized to allow
them to be used with simple data structures and with an indexible
data structure like packcontainer.
Wavelet packet tree code
The packtree and packfreq classes (derived from
the common base class packtree_base) build wavelet packet
trees. In the case of the packfreq class, the nodes at a given
level (a so called level basis) of the tree are ordered by frequency,
(moving from left to right).
The wavelet packet algorithm is independent of the actual wavelet
algorithm that is used while building the tree. The following
wavelet algorithms are included:
-
Daubechies D4. The forward and inverse Daubechies D4 algorith is
supported for building a standard wavelet packet tree (which can
be used in best basis calculation). Only the forward wavelet
step functions are provided for calculating a frequency ordered
wavelet packet tree.
-
Haar
-
Lifting scheme version of the Haar wavelet (see
Basic
Lifting Scheme Wavelets).
-
Haar "classic". This is the most common version of the Haar
wavelet.
-
Haar "classic" for frequency ordering. This includes the forward
and inverse wavelets steps for calculating a frequency ordered
wavelet packet tree and for calculating the inverse.
-
Line wavelet. This is a lifting scheme wavelet that uses linear
interpolation.
The lossless wavelet packet tree code uses integer to integer versions
of the Haar and line wavelet algorithms. Another algorithm, sometimes
called the TS transform, which is popular in image processing is also
included.
Once the wavelet packet tree is constructed a "best basis" algorithm
can be applied. The "best basis" is the minimal representation of
the data relative to a particular cost function. This code
provides cost functions derived from a common base class,
costbase.
The packcontainer class allows a sequence of N
numbers to be treated as both a contiguous array or as a left
half and a right half (each consisting of N/2 data elements).
The wavelet packet and modified wavelet packet algorithms were
developed for floating point (e.g., double) arrays (or containers with
double elements). The lossless wavelet packet compression code uses
integer to integer transforms. In writing this code I have attempted
to share as much as possible with the original floating point wavelet
packet code. However, in some cases the lossless compression code
consists of integer specific versions of the original wavelet packet
code. For example, costbase_int is the integer specific cost
function base class. The packtree_int and
invpacktree_int classes provide integer specific wavelet
packet tree support.
In most cases memory is allocated from a memory pool. This
approach to dynamic allocation allows memory to be allocated
without worrying about deallocation. At a point in the program
when the memory is no longer needed, it is all deallocated
at once by deallocating the memory pool. The only exception to
pool allocation is the GrowableArray class, which uses the standard
new and delete operators.
Copyright and Use
You may use this source code without limitation and without
fee as long as you include:
This software was written and is copyrighted by Ian Kaplan, Bear
Products International, www.bearcave.com, 2002.
This software is provided "as is", without any warranty or
claim as to its usefulness. Anyone who uses this source code
uses it at their own risk. Nor is any support provided by
Ian Kaplan and Bear Products International.
Please send any bug fixes or suggested source changes to:
iank@bearcave.com
*/
packet/Makefile 0100644 0000764 0000764 00000006154 07525020700 012427 0 ustar iank iank
#
# Makefile for Wavelet Packet code
#
# This makefile assumes that the Microsoft environment has
# been enhanced with the cygwin UNIX utilities for Win32.
# These are available from Red Hat (which purchased Cygnus).
#
# Command:
# nmake -f Makefile
#
DEBUG = -Zi
BROWSE = -FR
CFLAGS = $(BROWSE) $(DEBUG) -DWIN32 -Tp
OBJ = obj
EXE = .exe
LOCALINC = include
DOXYPATH = e:\doxygen\bin
INCLUDE = -I$(LOCALINC) -I "d:\Program Files\Microsoft Visual Studio\Vc98\Include"
PACKET_OBJ = blockpool.$(OBJ) packtree_base.$(OBJ) packtree.$(OBJ) packfreq.$(OBJ) costshannon.$(OBJ) invpacktree.$(OBJ)
all: packtest freqtest
packtest: packtest$(EXE)
freqtest: freqtest$(EXE)
packtest$(EXE): $(PACKET_OBJ) packtest.$(OBJ)
$(CC) $(PACKET_OBJ) packtest.$(OBJ) $(DEBUG) -o packtest$(EXE)
freqtest$(EXE): $(PACKET_OBJ) freqtest.$(OBJ)
$(CC) $(PACKET_OBJ) freqtest.$(OBJ) $(DEBUG) -o freqtest$(EXE)
#
# clean-up for Microsoft object and browser files and emacs temps
#
clean:
rm -f *.obj *.pdb *.sbr *.ilk *.exe
rm -f include/*~
rm -f src/*~
doxygen:
$(DOXYPATH)\doxygen doxygenDocConfig
include\costbase.h: include\packnode.h
include\costshannon.h: include\costbase.h include\packnode.h
include\costthresh.h: include\costbase.h include\packnode.h
include\haar.h: include\liftbase.h
include\haar_classic.h: include\liftbase.h
include\line.h: include\liftbase.h
include\packdata.h: include\blockpool.h
include\invpacktree.h: include\packdata.h
include\packnode.h: include\packdata.h
include\packdata_list.h: include\fifo_list.h include\packdata.h
include\packcontainer.h: include\packnode.h
include\queue.h: include\fifo_list.h include\packnode.h
include\packtree_base.h: include\packnode.h include\liftbase.h
include\packtree.h: include\packtree_base.h include\packcontainer.h include\liftbase.h
include\packfreq.h: include\packtree_base.h include\liftbase.h
include\invpacktree.h: include\liftbase.h include\list.h include\packcontainer.h include\packdata.h include\packdata_list.h
packtree_base.$(OBJ): src\$*.cpp include\$*.h
$(CC) -c $(INCLUDE) $(CFLAGS) src\$*.cpp
packtree.$(OBJ): src\$*.cpp include\$*.h include\queue.h include\packcontainer.h include\blockpool.h
$(CC) -c $(INCLUDE) $(CFLAGS) src\$*.cpp
packfreq.$(OBJ): src\$*.cpp include\$*.h
$(CC) -c $(INCLUDE) $(CFLAGS) src\$*.cpp
costshannon.$(OBJ): src\$*.cpp include\$*.h
$(CC) -c $(INCLUDE) $(CFLAGS) src\$*.cpp
invpacktree.$(OBJ): src\$*.cpp include\$*.h include\blockpool.h
$(CC) -c $(INCLUDE) $(CFLAGS) src\$*.cpp
blockpool.$(OBJ): src\$*.cpp include\$*.h
$(CC) -c $(INCLUDE) $(CFLAGS) src\$*.cpp
local_new.$(OBJ): src\$*.cpp
$(CC) -c $(INCLUDE) $(CFLAGS) src\$*.cpp
packtest.$(OBJ): src\$*.cpp include\haar.h include\haar_classic.h include\packnode.h include\packcontainer.h include\packtree.h include\invpacktree.h include\costshannon.h include\costthresh.h
$(CC) -c $(INCLUDE) $(CFLAGS) src\$*.cpp
freqtest.$(OBJ): src\$*.cpp include\haar_classicFreq.h include\packcontainer.h
$(CC) -c $(INCLUDE) $(CFLAGS) src\$*.cpp
packet/Makefile.unix 0100644 0000764 0000764 00000005201 10176102663 013406 0 ustar iank iank
#
# Makefile for Wavelet Packet code on UNIX (Live Free or Die!)
#
# This make file assumes that the GNU g++ compiler is available
#
# To run this make file (on UNIX) enter
#
# make -f Makefile.unix
#
# Most recently tested using
#
# g++ (GCC) 3.2.3 20030502 (Red Hat Linux 3.2.3-42)
#
DEBUG = -g
BROWSE =
GNU_FLAGS = -Wno-deprecated
CFLAGS = $(BROWSE) $(DEBUG) -D_UNIX_ $(GNU_FLAGS)
OBJ = o
#
# Specify the UNIX C++ compiler here
#
CC = g++
# CC = CC
# CC = /usr/local/forte/SUNWspro/bin/CC
#
LOCALINC = include
INCLUDE = -I$(LOCALINC)
PACKET_OBJ = blockpool.$(OBJ) packtree_base.$(OBJ) packtree.$(OBJ) packfreq.$(OBJ) costshannon.$(OBJ) invpacktree.$(OBJ) packtest.$(OBJ)
all: $(PACKET_OBJ)
$(CC) $(PACKET_OBJ) $(DEBUG) -o packet
#
# clean-up for Microsoft object and browser files and emacs temps
#
# It's totally wierd, clean stopped working...
clean:
# rm *.obj *.pdb *.sbr *.exe
# rm -f include/*~
# rm -f src/*~
include/costbase.h: include/packnode.h
include/costshannon.h: include/costbase.h include/packnode.h
include/costthresh.h: include/costbase.h include/packnode.h
include/haar.h: include/liftbase.h
include/haar_classic.h: include/liftbase.h
include/line.h: include/liftbase.h
include/packdata.h: include/blockpool.h
include/invpacktree.h: include/packdata.h
include/packnode.h: include/packdata.h
include/packdata_list.h: include/fifo_list.h include/packdata.h
include/packcontainer.h: include/packnode.h
include/queue.h: include/fifo_list.h include/packnode.h
include/packtree.h: include/packnode.h include/packcontainer.h include/liftbase.h
include/invpacktree.h: include/liftbase.h include/list.h include/packcontainer.h include/packdata.h include/packdata_list.h
packtree_base.$(OBJ): include/queue.h include/packnode.h include/packcontainer.h include/packtree_base.h
$(CC) -c $(INCLUDE) $(CFLAGS) src/$*.cpp
packtree.$(OBJ): include/queue.h include/packcontainer.h include/blockpool.h
$(CC) -c $(INCLUDE) $(CFLAGS) src/$*.cpp
packfreq.$(OBJ): include/packfreq.h
$(CC) -c $(INCLUDE) $(CFLAGS) src/$*.cpp
costshannon.$(OBJ):
$(CC) -c $(INCLUDE) $(CFLAGS) src/$*.cpp
packdata_list.$(OBJ):
$(CC) -c $(INCLUDE) $(CFLAGS) src/$*.cpp
invpacktree.$(OBJ): include/blockpool.h
$(CC) -c $(INCLUDE) $(CFLAGS) src/$*.cpp
blockpool.$(OBJ):
$(CC) -c $(INCLUDE) $(CFLAGS) src/$*.cpp
local_new.$(OBJ):
$(CC) -c $(INCLUDE) $(CFLAGS) src/$*.cpp
packtest.$(OBJ): include/haar.h include/haar_classic.h include/packnode.h include/packcontainer.h include/packtree.h include/invpacktree.h include/costshannon.h include/costthresh.h
$(CC) -c $(INCLUDE) $(CFLAGS) src/$*.cpp
%.o: src/%.cpp include/%.h
packet/src/ 0040755 0000764 0000764 00000000000 10176075422 011562 5 ustar iank iank packet/src/blockpool.cpp 0100644 0000764 0000764 00000015550 10176075422 014255 0 ustar iank iank
/** \file
This file contains the class functions for a "pool" based memory
allocator. These functions manage the block chain that makes up the
memory pool and allocates memory from the individual blocks. The
entire pool is freed at once.
The pool allocator was originally written as a low level allocation
package. The original package queries the system to find the page
size and allocation block size. This version has been simplified,
but some of the low level nature remains.
The documentation in this file is formatted for doxygen
(see www.doxyeng.org).
Copyright and Use
You may use this source code without limitation and without
fee as long as you include:
This software was written and is copyrighted by Ian Kaplan, Bear
Products International, www.bearcave.com, 2002.
This software is provided "as is", without any warranty or
claim as to its usefulness. Anyone who uses this source code
uses it at their own risk. Nor is any support provided by
Ian Kaplan and Bear Products International.
Please send any bug fixes or suggested source changes to:
iank@bearcave.com
@author Ian Kaplan
*/
#include
#include
#include
#include "blockpool.h"
unsigned int block_pool::alloc_gran = (unsigned int)block_pool::page_size;
block_pool::block_chain *block_pool::current_block = 0;
block_pool::block_chain *block_pool::block_list_start = 0;
/**
block_pool::init_pool
This function is automatically called to initialize the block_pool
object when the "current_block" pointer is null.
*/
void block_pool::init_pool( void )
{
block_chain *new_link;
new_link = new_block( alloc_gran );
block_list_start = new_link;
current_block = new_link;
} /* init_pool */
/**
new_block
The new_block function is the "root" memory allocator for the
block_pool object. The amount of memory allocated is rounded up to
the next "block_size" boundary. Both the block_chain structure and
the allocatible memory are allocated from a single block that is a
multiple of the page size. This should avoid fragmentation in the
system memory allocator.
The "page" referenced here is a virtual memory page. This version
of the code does not actually check the system page size, but assumes
a 4Kb page.
*/
block_pool::block_chain *block_pool::new_block( unsigned int block_size )
{
const unsigned int max_block_size = max_block_multiple * page_size;
block_chain *new_link = 0;
unsigned int alloc_amt, total_alloc;
// add in the memory needed for the block_chain structure
total_alloc = block_size + sizeof(block_chain);
if (total_alloc < alloc_gran)
alloc_amt = alloc_gran;
else {
// its larger than the minimum allocation granularity, so round
// up the the nearest page.
alloc_amt = ((total_alloc + (page_size-1))/page_size) * page_size;
}
if (alloc_amt <= max_block_size) {
/* Allocate memory for both the block_chain structure and the memory
block */
new_link = (block_chain *)MemAlloc( alloc_amt );
// The new memory block starts after the block_chain structure
Chain_block(new_link) = (void *)(((unsigned int)new_link) + sizeof(block_chain));
assert( alloc_amt >= block_size );
Chain_bytes_used(new_link) = 0;
Chain_block_size(new_link) = alloc_amt - sizeof(block_chain);
Chain_next(new_link) = 0;
}
else {
printf("block_pool::new_block: allocation request too large\n");
}
return new_link;
} // block_chain::new_block
/**
block_pool::add_block Add a new memory block to the memory pool.
This function is called when the amount of memory requested by
pool_alloc will not fit in the current block.
*/
void *block_pool::add_block( unsigned int block_size )
{
block_chain *block = 0;
block_chain *last_block;
last_block = current_block;
block = new_block( block_size );
Chain_next(current_block) = block;
current_block = block;
return (void *)block;
} // block_chain::add_block
/**
pool_alloc
Allocate memory from the memory pool. The pool_alloc and free_pool
functions do memory allocation and deallocation.
This function is called to allocate memory from the memory pool.
If there is enough free memory in the current block to satisify the
memory request, memory is allocated from the current block and the
amount of free memory is updated. If the current block does not
have enough memory, add_block is called to allocate a new memory
block which will be large enough.
*/
void *block_pool::pool_alloc( unsigned int num_bytes )
{
const unsigned int align = sizeof( int );
void *addr = 0;
unsigned int amt_free;
/* the number of bytes allocated must be a multiple of the align
size */
num_bytes = ((num_bytes + (align-1))/align) * align;
if (current_block == 0) {
init_pool();
}
amt_free = Chain_block_size(current_block) - Chain_bytes_used(current_block);
if (num_bytes > amt_free) {
if (add_block( num_bytes ) != 0) {
amt_free = Chain_block_size(current_block);
}
}
if (amt_free >= num_bytes) {
addr = (void *)((unsigned int)Chain_block(current_block) + Chain_bytes_used(current_block));
Chain_bytes_used(current_block) += num_bytes;
}
else {
printf("block_pool::block_alloc: allocation error\n");
exit(1);
}
return addr;
} // block_pool::pool_alloc
/**
block_pool::free_pool
Walk through the block chain and deallocate the blocks.
Note that the block chain structures and the allocatible
memory is contained within a single allocated block.
The block_chain structure is at the start of this block
so passing its address to the memory deallocation function
deallocates both the block chain structure and the allocatible
memory
*/
void block_pool::free_pool(void)
{
block_chain *tmp;
while (block_list_start != 0) {
tmp = block_list_start;
block_list_start = Chain_next(block_list_start);
MemFree( (void *)tmp );
}
} // free_pool
/**
print_block_pool_info
Print information about the block pool
*/
void block_pool::print_block_pool_info( FILE *fp /*= stdout */)
{
int total_allocated = 0;
int total_unused = 0;
block_chain *ptr = block_list_start;
fprintf(fp, "Minimum memory allocation size: %d\n", alloc_gran );
fprintf(fp, "Page size: %d\n", (unsigned int)page_size );
fprintf(fp, "[block size, bytes_used]\n");
while (ptr != 0) {
fprintf(fp, "[%4d, %4d]", Chain_block_size(ptr), Chain_bytes_used(ptr));
total_allocated += Chain_bytes_used(ptr);
total_unused += (Chain_block_size(ptr) - Chain_bytes_used(ptr));
if (Chain_next(ptr) != 0) {
fprintf(fp, ", ");
}
else {
fprintf(fp, "\n");
}
ptr = Chain_next(ptr);
} // while
fprintf(fp, "Total allocated = %5d, total unused = %3d\n", total_allocated,
total_unused );
}
packet/src/costshannon.cpp 0100644 0000764 0000764 00000001554 07515164511 014626 0 ustar iank iank
#include
#include
#include
#include "costshannon.h"
/**
An implementation of a modified version of the Shannon entropy
function as a wavelet packet cost function. This is described
in section 8.3.2 of Ripples in Mathematics by
Jensen and la Cour-Harbo.
The log function here is the natural log (sometimes denoted as
ln()). Note that the result of the entropy function is always
negative.
*/
double costshannon::costCalc( packnode *node )
{
assert( node != 0 );
size_t len = node->length();
const double *a = node->getData();
double sum = 0.0;
for (int i = 0; i < len; i++) {
double val = 0.0;
if (a[i] != 0.0) {
double square = a[i] * a[i];
val = square * log( square );
}
sum = sum + val;
}
return -sum;
} // costshannon
packet/src/freqtest.cpp 0100644 0000764 0000764 00000013146 07525022502 014120 0 ustar iank iank
/** \file
This file contains code to test wavelet frequency analysis via
wavelet packets. See
http://www.bearcave.com/misl/misl_tech/wavelets/packfreq/index.html
The documentation in this file is formatted for doxygen
(see www.doxygen.org).
Copyright and Use
You may use this source code without limitation and without
fee as long as you include:
This software was written and is copyrighted by Ian Kaplan, Bear
Products International, www.bearcave.com, 2002.
This software is provided "as is", without any warranty or
claim as to its usefulness. Anyone who uses this source code
uses it at their own risk. Nor is any support provided by
Ian Kaplan and Bear Products International.
Please send any bug fixes or suggested source changes to:
iank@bearcave.com
@author Ian Kaplan
*/
#include
#include
#include "packcontainer.h"
#include "haar_classicFreq.h"
#include "daub.h"
#include "packfreq.h"
/**
Some small data sets
*/
/**/
double data[] = { 32.0, 10.0, 20.0, 38.0, 37.0, 28.0, 38.0, 34.0,
18.0, 24.0, 18.0, 9.0, 23.0, 24.0, 28.0, 34.0 };
/**/
/*
double data[] = { 56.0, 40.0, 8.0, 24.0, 48.0, 48.0, 40.0, 16.0 };
*/
/** \function
Generate a signal composed of a sum of sine waves.
The sine waves summed are in decreasing magnitude and
decreasing frequency.
*/
void gen_freqMix( double *vecX, double *vecY, size_t N)
{
const double PI = 3.1415926535897932384626433832795;
const double range = 2 * PI;
const double incr = range / (double)N;
double point = 0.0;
int i;
for (i = 0; i < N; i++) {
vecX[i] = point;
vecY[i] = 4 * sin( 64 * point ) +
2 * sin( 32 * point ) +
1 * sin( 16 * point ) +
0.5 * sin( 8 * point );
// printf("%7.4g, %7.4g\n", vecX[i], vecY[i] );
point = point + incr;
}
} // gen_freqMix
/** \function
Generate a signal composed of one or more sine waves.
*/
void gen_sinCombo( double *vecX, double *vecY, size_t N )
{
const double PI = 3.1415926535897932384626433832795;
const double range = 8 * PI;
const double incr = range / (double)N;
double point = 0.0;
int stepCnt = 0;
int i;
for (i = 0; i < N; i++) {
vecX[i] = point;
// vecY[i] = sin( 64 * point ) + sin( 32 * point ) + sin( 16 * point );
// vecY[i] = sin( 16 * PI * point ) + sin( 4 * PI * point );
vecY[i] = sin( 4 * PI * point );
// printf("x[%2d] = %7.4g, y[%2d] = %7.4g\n", i, vecX[i], i, vecY[i] );
// printf("%7.4g, %7.4g\n", vecX[i], vecY[i] );
// printf("%d, %7.4g\n", i, vecY[i] );
point = point + incr;
}
}
/** \function
Generate a signal that increases in frequency by steps.
*/
void gen_steps( double *vecX, double *vecY, size_t N, size_t steps )
{
const double PI = 3.1415926535897932384626433832795;
const double range = 32 * PI;
const double incr = range / (double)N;
double point = 0.0;
double mult = 1;
const size_t stepWidth = N / steps;
int stepCnt = 0;
int i;
for (i = 0; i < N; i++) {
vecX[i] = point;
vecY[i] = sin( mult * point );
// printf("x[%2d] = %7.4g, y[%2d] = %7.4g\n", i, vecX[i], i, vecY[i] );
// printf("%7.4g, %7.4g\n", vecX[i], vecY[i] );
// printf("%d, %7.4g\n", i, vecY[i] );
point = point + incr;
stepCnt++;
if (stepCnt == stepWidth) {
mult = mult + (PI/2.0);
stepCnt = 0;
}
}
}
/** \function
Generate a linear chirp signal
*/
void gen_chirp( double *vecX, double *vecY, size_t N )
{
const double PI = 3.1415926535897932384626433832795;
const double range = 2;
// const double range = 16 * PI;
const double incr = range / (double)N;
double point = 0.0;
int i;
for (i = 0; i < N; i++) {
vecX[i] = point;
vecY[i] = sin( 128 * PI * point * point );
// printf("x[%2d] = %7.4g, y[%2d] = %7.4g\n", i, vecX[i], i, vecY[i] );
// printf("%7.4g, %7.4g\n", vecX[i], vecY[i] );
// printf("%d, %7.4g\n", i, vecY[i] );
point = point + incr;
}
}
/** \function
Print a vector of doubles whose length is len.
*/
void prCoords( double *vecX, double *vecY, size_t len )
{
for (int i = 0; i < len; i++) {
printf("%7.4f %7.4f\n", vecX[i], vecY[i] );
}
}
/** \function
Print a vector of doubles
*/
void prVec( double *vec, size_t len )
{
for (int i = 0; i < len; i++) {
printf("%4d %7.4f\n", i, vec[i] );
}
} // prVec
/** \function
The entry point for code to test the wavelet packet transform.
The code in main provides a simple example of how to call the
wavelet packet transform code.
*/
int
main()
{
const size_t N = 1024;
// const size_t N = sizeof( data ) / sizeof( double );
double vecX[N], vecY[N];
// gen_chirp( vecX, vecY, N );
// gen_steps( vecX, vecY, N, 8 );
gen_freqMix( vecX, vecY, N );
// gen_sinCombo( vecX, vecY, N );
// prVec( vecY, N );
// prCoords( vecX, vecY, N );
// The "Haar" classic transform
haar_classicFreq h;
// Daubechies h;
// calculate the wavelet packet tree, using the wavelet transform h
packfreq tree( vecY, N, &h );
// packfreq tree( data, N, &h );
// tree.pr();
// printf("\n");
tree.getLevel( 5 );
tree.plotMat(N);
// free the memory pool
block_pool mem_pool;
mem_pool.free_pool();
return 0;
}
packet/src/invpacktree.cpp 0100644 0000764 0000764 00000015247 07515164710 014611 0 ustar iank iank
#include
#include
#include "blockpool.h"
#include "invpacktree.h"
/** \file
For additional information on this implementation of the inverse
wavelet packet transform see
http://www.bearcave.com/misl/misl_tech/wavelets/packet/index.html
The documentation in this file is formatted for doxygen
(see www.doxyeng.org).
Copyright and Use
You may use this source code without limitation and without
fee as long as you include:
This software was written and is copyrighted by Ian Kaplan, Bear
Products International, www.bearcave.com, 2002.
This software is provided "as is", without any warranty or
claim as to its usefulness. Anyone who uses this source code
uses it at their own risk. Nor is any support provided by
Ian Kaplan and Bear Products International.
Please send any bug fixes or suggested source changes to:
iank@bearcave.com
@author Ian Kaplan
*/
/**
Allocate a new level in the wavelet packet tree that is being
constructed. The new "level" is built from a packcontainer object.
The left child will be the elem argument. The length of the data in
the new level will be twice that of lower level (of which elem is
half).
*/
void invpacktree::new_level( packdata *elem )
{
size_t half = elem->length();
size_t n = half * 2;
packcontainer *container = new packcontainer( n );
container->lhsData( (double *)elem->getData() );
stack.add( container );
} // new_level
/**
At this point the Top Of Stack (TOS) packcontainer object should
have both a right and a left array. Calculate an inverse wavelet
transform step on the packcontainer object. The packcontainer
object allows the right and left hand side arrays to be treated as
one array
If the current top of stack is twice the size of the inverse wavelet
transform step result, the result becomes the right hand size of
top of stack packcontainer and reduce is called recursively.
If the TOS is empty or it is not twice the size of the inverse
transform result, a new packcontainer will be pushed on the stack.
The left hand size will be the transform result.
*/
void invpacktree::reduce()
{
LIST::handle h;
h = stack.first();
packcontainer *tos = stack.get_item( h );
assert( tos->lhsData() != 0 && tos->rhsData() != 0 );
/** Remove the linked list element that used to contain tos
(e.g., pop the linked list element off). Note that this
leaves the tos object unchanged (e.g., it does not delete
it).
*/
stack.remove();
size_t n = tos->length();
// calculate the inverse wavelet transform step
waveObj->inverseStep( (*tos), n );
// copy the result of the inverse wavelet transform
// into a new data array.
block_pool mem_pool;
double *vec = (double *)mem_pool.pool_alloc( n * sizeof( double ) );
for (int i = 0; i < n; i++) {
vec[i] = (*tos)[i];
}
if (stack.first() != 0) {
h = stack.first();
packcontainer *tos = stack.get_item( h );
if (tos->length() == n*2) {
tos->rhsData( vec );
reduce();
}
else {
assert( tos->length() > n*2 );
packcontainer *container = new packcontainer( n*2 );
container->lhsData( vec );
stack.add( container );
} // else
}
else {
// the stack is empty
packcontainer *container = new packcontainer( n*2 );
container->lhsData( vec );
stack.add( container );
}
} // reduce
/**
Add a packdata element to the inverse wavelet packet transform
calculation.
A packdata element contains the wavelet packet data items from
one level of the wavelet packet transform tree.
The inverse wavelet packet transform calculation uses a stack. Each
level in the stack consists of a packcontainer object. This object
consists of left and right hand size arrays. When both these arrays
are present, they can be treated as one contiguous array (courtesy
of C++ operator overloading). The size of the packcontainer object
is twice the size of its left and right hand size arrays.
If the stack is empty, a new level is created. A packetcontainer
object can be viewed as a binary tree element. The left hand
size is filled in first.
If the stack is not empty and the packcontainer object on the
top of stack (TOS) is twice the size of the elem argument,
then the array contained in the elem argument is added to
the TOS element and reduce is called to calculate a step
of the inverse wavelet transform.
If the TOS element is greater than twice the size of elem
then a new level is added.
*/
void invpacktree::add_elem( packdata *elem )
{
assert( elem != 0 );
if (stack.first() == 0) {
new_level( elem );
}
else {
size_t half = elem->length();
size_t n = half * 2;
LIST::handle h;
h = stack.first();
packcontainer *tos = stack.get_item( h );
if (tos->length() == n) {
assert( tos->rhsData() == 0);
tos->rhsData( (double *)elem->getData() );
reduce();
}
else if (tos->length() > n) {
new_level( elem );
}
else {
printf("add_elem: the size of the TOS elem is wrong\n");
}
} // else
} // add_elem
/**
This constructor calculates the inverse wavelet packet transform.
The constructor is passed a packdata_list object, which contains
the "best basis" result of the wavelet packet transform. The
liftbase template argument is a pointer to a wavelet transform
function. This wavelet transform must be the same function that
was used to calculate the packet transform.
*/
invpacktree::invpacktree( packdata_list &list,
liftbase *w )
{
data = 0;
N = 0;
waveObj = w;
// Traverse the "best basis" list and calculate the inverse
// wavelet packet transform.
packdata_list::handle h;
for (h = list.first(); h != 0; h = list.next( h )) {
packdata *elem = list.get_item( h );
add_elem( elem );
} // for
LIST::handle tosHandle;
tosHandle = stack.first();
packcontainer *tos = stack.get_item( tosHandle );
if (tos != 0) {
size_t len = tos->length();
N = len/2;
data = tos->lhsData();
stack.remove();
}
} // invpacktree
/**
print the result of the inverse wavelet packet transform
*/
void invpacktree::pr()
{
if (data != 0) {
for (int i = 0; i < N; i++) {
printf("%7.4f ", data[i] );
}
printf("\n");
}
} // pr
packet/src/local_new.cpp 0100644 0000764 0000764 00000001674 07466637225 014253 0 ustar iank iank
#include
#include
/** \file
This file contains overrides for the global new and delete
operators. These exist to assure that only the pool allocation
functions are called. If you see a message from one of these
functions then something other than pool allocation is taking
place.
This file is only used for testing. If you don't want to do this
check you can remove this file from the software build.
*/
void *operator new( size_t num_bytes )
{
printf("global operator new\n");
void *rtn = malloc( num_bytes );
return rtn;
} // new
void *operator new[]( size_t num_bytes )
{
printf("global operator new []\n");
void *rtn = malloc( num_bytes );
return rtn;
}
void operator delete( void *addr )
{
printf("global operator delete\n");
free( addr );
}
void operator delete[](void *addr )
{
printf("global operator delete []\n");
free( addr );
}
packet/src/packfreq.cpp 0100644 0000764 0000764 00000021174 07515164006 014064 0 ustar iank iank
/** \file
The documentation in this file is formatted for doxygen
(see www.doxyeng.org).
Copyright and Use
You may use this source code without limitation and without
fee as long as you include:
This software was written and is copyrighted by Ian Kaplan, Bear
Products International, www.bearcave.com, 2002.
This software is provided "as is", without any warranty or
claim as to its usefulness. Anyone who uses this source code
uses it at their own risk. Nor is any support provided by
Ian Kaplan and Bear Products International.
Please send any bug fixes or suggested source changes to:
iank@bearcave.com
@author Ian Kaplan
*/
#include
#include
#include
#include "packfreq.h"
/**
Construct a wavelet packet tree from a vector of double values.
The size of the vector, which must be a power of two, is passed in
N. A pointer to a wavelet Lifting Scheme object is passed in
w. The wavelet Lifting Scheme object is used to calculate
the wavelet transform step which is applied at each level (where
level > 0) of the wavelet packet tree.
The first level (level 0) of the wavelet packet tree contains
the original data set.
The "newLevel" function is defined in the base class. This
function is used to calculate both the standard and the
modified wavelet packet trees. The modified wavelet packet
tree is calculated by this class and is used for Wavelet
frequency analysis.
The wavelet transform step uses two filters H and G to
calculate the low pass (scaling function) and the high pass (wavelet
function). In the standard transform the low pass result is always
placed in the lower half of the array and the high pass result is
placed in the upper half of the array. Each of these two half
arrays form the input for the next step. In the modified transform
calculated by this class, the right child uses the standard
algorithm where H and G results are in the upper and
lower halves of the result array. However, the left child inverts
the locations of the H and G results, placing the
result of the G filter in the lower half and the result
of the H filter in the upper half.
The two boolean flags that are passed to the newLevel function
determine whether a frequency ordered wavelet packet tree
is calculated and whether the location of the filter results
is inverted.
\arg vec An array of double values on which the wavelet packet
transform is calculated.
\arg N The number of elements in the input array
\arg w A pointer to the the wavelet transform object to use
in calculating the wavelet packet transform.
*/
packfreq::packfreq( const double *vec,
const size_t N,
liftbase *w )
{
waveObj = w;
block_pool mem_pool;
double *vecCopy = (double *)mem_pool.pool_alloc( N * sizeof( double ) );
for (int i = 0; i < N; i++) {
vecCopy[i] = vec[i];
}
root = new packnode( vecCopy, N, packnode::OriginalData );
root->mark( true );
//
// The first level uses the standard wavelet calculation, so
// reverse = false
// freqCalc, reverse
newLevel( root, true, false );
} // packfreq
/**
Traverse the tree, from left to right, and add
a node at level to the level basis matrix.
Levels are numbered from 0 at the root.
*/
void packfreq::findLevel( packnode* top,
size_t cur_level,
const size_t level )
{
if (top != 0) {
if (cur_level == level) {
mat.append(top);
}
else {
findLevel( top->lhsChild(), cur_level+1, level );
findLevel( top->rhsChild(), cur_level+1, level );
}
}
} // findLevel
/**
Build a level basis matrix
The "level basis" is a horizontal slice through the wavelet packet
tree. Here the levels are numbered from zero at root (which
contains the original data). The most useful basis for frequency
analysis is a square matrix. If the original data set consists of
1024 elements, level 5 of the tree consists of 32 tree nodes
containing 32 elements each.
number of elements
in a tree node level number of nodes
1024 0 1
512 1 2
256 2 4
128 3 8
64 4 16
32 5 32
The level basis is built from left to right, in the horizontal slice
through the tree. The left most node contains the lowest frequency
band.
*/
void packfreq::getLevel( const size_t level )
{
findLevel( root, 0, level );
} // getLevel
/**
Print out the level basis matrix so that it can be plotted as a
three dimensional surface.
The level basis for a given level is constructed by the getLevel()
function.
In the discussion of this function below, we will consider the
case where the modified wavelet packet tree was built from a
data set of 1024 elements and the level basis was at level
5, resulting in a matrix of 32 nodes, each of which has
32 values.
mat[0] is the left most node (and lowest frequency band), mat[31] is
the right most (and highest frequency band). The values
(*mat[0])[0] ... (*mat[0])[31] represent frequency values at
particular time intervals.
If we use two variables, x and y, to index the matrix then
(*mat[y])[x] defines a time/frequency plane, where the x-index
is time and the y-index is frequency. The value at (*mat[y])[x]
is the magnitude a particular time/frequency point.
In plotting a surface from this data, time or frequency are plotted
on the x or y axis. The magnitude at (*mat[y])[x] is plotted on the
Z-axis. Following Ripples in Mathematics by Jensen and
la-Cour Harbo (Springer Verlag, 2001), the calculation for the
z-axis value is:
m = (*mat)[y][x]
z_val = log( 1 + val^2 )
Where log is the natural log.
This function has been used to generate data for gnuPlot surface
plots.
Whether frequency is plotted on the x-axis or the y-axis depends
on the nature of the data. If a signal line sin(x) is analyzed,
the frequency is constant. To avoid the ridge produced by the
sin(x) frequency obscuring the surface, it is better to plot
time on the y-axis and frequency on the x-axis.
In the case of a constantly changing frequency (like the linear
chirp), time is plotted on the x-axis and frequency is plotted
on the y-axis.
Another wrinkle involves the values plotted on the x and y axis.
The simplest plot of time (say on the x-axis) would consist of the
32 time regions (e.g., (*mat[y][x])). Similarly, frequency would
consist of 32 frequency regions (on the y-axis). This would
yield a surface numbered from 0..31 on both the x and y axis.
The axis labels can be scaled to represent the frequency and
time ranges.
This function was used to generate different plots, so various parts
are commented out, making it a bit of a hack.
*/
void packfreq::plotMat(const size_t N)
{
size_t num_y = mat.length();
const double incr = (double)N / (double)num_y;
if (num_y > 0) {
size_t num_x = mat[0]->length();
double freq_start = 0.0;
for (size_t y = 0; y < num_y; y++) {
double time_start = 0.0;
for (size_t x = 0; x < num_x; x++) {
double val = (*mat[y])[ x ];
// plot time on the x, frequency on y
// printf(" %d %d %7.4f\n", x, y, log(1+(val*val)) );
// plot frequency on x, time on y
printf(" %d %d %7.4f\n", y, x, log(1+(val*val)) );
// plot actual frequency and time values,
// with frequency on x and time on y
// printf(" %7.4f %7.4f %7.4f\n",
// freq_start, time_start, log(1+(val*val)) );
time_start = time_start + incr;
}
freq_start = freq_start + incr;
printf("\n");
}
}
} // plotMat
/**
Print the contents of the level basis matrix
*/
void packfreq::prMat()
{
int num_y = mat.length();
if (num_y > 0) {
size_t num_x = mat[0]->length();
for (int y = num_y-1; y >= 0; y--) {
for (size_t x = 0; x < num_x; x++) {
printf(" %7.4f ", (*mat[y])[ x ] );
}
printf("\n");
fflush(stdout);
}
}
} // prMat
packet/src/packtest.cpp 0100644 0000764 0000764 00000010652 07525021700 014077 0 ustar iank iank
/** \file
This file contains test code for wavelet packet tree construction
and for testing the best basis calculation.
The documentation in this file is formatted for doxygen
(see www.doxygen.org).
Copyright and Use
You may use this source code without limitation and without
fee as long as you include:
This software was written and is copyrighted by Ian Kaplan, Bear
Products International, www.bearcave.com, 2002.
This software is provided "as is", without any warranty or
claim as to its usefulness. Anyone who uses this source code
uses it at their own risk. Nor is any support provided by
Ian Kaplan and Bear Products International.
Please send any bug fixes or suggested source changes to:
iank@bearcave.com
@author Ian Kaplan
*/
#include
#include
#include "haar.h"
#include "haar_classic.h"
#include "haar_classicFreq.h"
#include "packnode.h"
#include "packcontainer.h"
#include "packtree.h"
#include "invpacktree.h"
#include "costshannon.h"
#include "costthresh.h"
/** \file
Test code for the floating point form of the wavelet packet
algorithm.
*/
double data[] = { 32.0, 10.0, 20.0, 38.0, 37.0, 28.0, 38.0, 34.0,
18.0, 24.0, 18.0, 9.0, 23.0, 24.0, 28.0, 34.0 };
/**
double data[] = { 56.0, 40.0, 8.0, 24.0, 48.0, 48.0, 40.0, 16.0 };
*/
/** \function
Print a vector of doubles whose length is len.
*/
void prVec( double *vec, size_t len )
{
for (int i = 0; i < len; i++) {
printf("%7.4f ", vec[i] );
}
printf("\n");
}
/**
Test that the wavelet transform is invertable
*/
void testWaveletTrans(const double *data, const size_t len )
{
haar haarVec;
block_pool mem_pool;
size_t num_bytes = len * sizeof(double);
double *vec = (double *)mem_pool.pool_alloc( num_bytes );
for (int i = 0; i < len; i++) {
vec[i] = data[i];
}
printf("Before forward trans:\n");
prVec( vec, len );
haarVec.forwardTrans( vec, len );
printf("After forward trans:\n");
prVec( vec, len );
haarVec.inverseTrans( vec, len );
printf("After inverse trans:\n");
prVec( vec, len );
printf("\n");
} // testWaveletTrans
/**
The entry point for code to test the wavelet packet transform.
The code in main provides a simple example of how to call the
wavelet packet transform code.
*/
int
main()
{
size_t len = sizeof(data)/sizeof(double);
// testWaveletTrans( data, len );
// The "Haar" classic transform
haar_classic h;
// calculate the wavelet packet tree, using the wavelet transform h
packtree tree( data, len, &h );
// print the wavelet transform tree (breadth first)
tree.pr();
// get the root of the wavelet packet transform tree
packnode *treeRoot = tree.getRoot();
// assign the Shannon entropy cost function to the tree
costshannon cost( treeRoot );
// Calculate a simple threshold cost function on the
// wavelet packet transform tree
// costthresh thresh( treeRoot, 1.0 );
printf("\n");
// Print the wavelet packet transform tree showing the cost
// function result.
tree.prCost();
// Calculate the "best basis" function from the tree
tree.bestBasis();
printf("\n");
// Print the wavelet packet tree showing the nodes selected
// as part of the "best basis" set.
tree.prBestBasis();
// Check that the best basis function succeeded. That is,
// that the best basis function does not include the
// original data.
if (tree.bestBasisOK()) {
printf("Best basis calculation succeeded\n");
}
else {
printf("Best basis calculation failed\n");
}
printf("\n");
// Get the best basis list. This will be a list of
// nodes consisting of the best basis set. This set is
// obtained by traversing the tree, top down, left to
// right.
packdata_list bestBasis = tree.getBestBasisList();
// Print the "best basis" set.
bestBasis.pr();
printf("\n");
// Calculate the inverse wavelet packet transform from the
// "best basis" list.
invpacktree invtree( bestBasis, &h );
printf("Inverse wavelet packet transform result:\n");
invtree.pr();
// free the memory pool
block_pool mem_pool;
mem_pool.free_pool();
return 0;
}
packet/src/packtree.cpp 0100644 0000764 0000764 00000017213 07524573350 014073 0 ustar iank iank
/** \file
The documentation in this file is formatted for doxygen
(see www.doxyeng.org).
Copyright and Use
You may use this source code without limitation and without
fee as long as you include:
This software was written and is copyrighted by Ian Kaplan, Bear
Products International, www.bearcave.com, 2002.
This software is provided "as is", without any warranty or
claim as to its usefulness. Anyone who uses this source code
uses it at their own risk. Nor is any support provided by
Ian Kaplan and Bear Products International.
Please send any bug fixes or suggested source changes to:
iank@bearcave.com
@author Ian Kaplan
*/
#include
#include
#include "packcontainer.h"
#include "packtree.h"
/**
Construct a wavelet packet tree from a vector of double
values. The size of the vector, which must be a power
of two, is passed in N. A wavelet Lifting Scheme object
is passed in w. This object is used to calculate
the wavelet transform step which is applied at each level
(where level > 0) of the wavelet packet tree.
The first level (level 0) of the wavelet packet tree contains
the original data set.
\arg vec An array of double values on which the wavelet packet
transform is calculated.
\arg N The number of elements in the input array
\arg w A pointer to the the wavelet transform object to use
in calculating the wavelet packet transform.
*/
packtree::packtree( const double *vec,
const size_t N,
liftbase *w )
{
waveObj = w;
block_pool mem_pool;
double *vecCopy = (double *)mem_pool.pool_alloc( N * sizeof( double ) );
for (int i = 0; i < N; i++) {
vecCopy[i] = vec[i];
}
root = new packnode( vecCopy, N, packnode::OriginalData );
root->mark( true );
newLevel( root, false, false );
} // packtree
/**
The best basis algorithm selects the nodes nearest the tree root for
the best basis set. These nodes are "marked" true with a boolean
flag.
The best basis algorithm outlined in Ripples in Mathematics
sets any marks in nodes that are below a marked node in the tree
(nodes which are in a subtree of a marked node) to false. However,
this is unnecessary when the best basis set is collected, since the
algorithm uses top down tree traversal and stops at marked node.
A problem does arise when the entire tree is printed to show the
nodes that are marked as part of the best basis set. In this case
the result may appear wrong, since child nodes of a best basis node
are marked. This function does a top down traversal setting the
"marks" on these child nodes to false.
*/
void packtree::cleanTree(packnode *top, bool removeMark )
{
if (top != 0) {
if (removeMark) {
if (top->mark()) {
top->mark( false );
}
}
else {
// if a mark is found, then set all the "marks" below that
// point to false (e.g., remove the marks).
if (top->mark()) {
removeMark = true;
}
}
cleanTree( top->lhsChild(), removeMark );
cleanTree( top->rhsChild(), removeMark );
}
} // cleanTree
/**
Print the wavelet packet tree cost values in breadth first
order.
*/
void packtree::prCost()
{
if (root != 0) {
breadthFirstPrint(printCost);
}
}
/**
Print the wavelet packet tree, showing the nodes
that have been selected by the "best basis" algorithm.
*/
void packtree::prBestBasis()
{
if (root != 0) {
cleanTree( root, false );
breadthFirstPrint(printBestBasis);
}
} // prBestBasis
/**
Walk the wavelet packet tree and apply the "best basis" algorithm
described in Chapter 8 of Ripples in Mathematics. Nodes that
are "marked" become part of the best basis, which is a minimal
representation of the data in terms of the cost function.
*/
double packtree::bestBasisWalk( packnode *top )
{
double cost = 0.0;
if (top != 0) {
packnode *lhs = top->lhsChild();
packnode *rhs = top->rhsChild();
if (lhs == 0 && rhs == 0) { // we've reached a leaf
cost = top->cost();
}
else if (lhs != 0 && rhs != 0) {
double lhsCost = bestBasisWalk( lhs );
double rhsCost = bestBasisWalk( rhs );
double v1 = top->cost();
double v2 = lhsCost + rhsCost;
if (v1 <= v2) {
top->mark( true );
lhs->mark( false );
rhs->mark( false );
}
else { // v1 > v2
top->cost( v2 );
}
cost = top->cost();
}
else {
// The tree does not seem to be a full binary tree
// Something has gone badly wrong.
assert( false );
}
}
return cost;
} // bestBasicWalk
/**
Calculate the wavelet packet "best basis"
*/
void packtree::bestBasis()
{
bestBasisWalk( root );
} // bestBasis
/**
Recursively traverse the wavelet packet tree and check that the best
basis result is correct. That is, that the best basis has been
calculated and that it does not include the orignal data set.
The best basis includes the original data set with the packnode
ofr the original data set is marked.
This algorithm makes use of two class global variables. There may
be a purely recursive, way to do this without these global class
variables, but these variables make the algorithm much easier.
The variables are initialized by the calling function
bestBasisOK().
*/
void packtree::checkBestBasis( packnode *top )
{
if (top != 0) {
if (top->mark()) {
foundBestBasisVal = true;
if (top->getKind() == packdata::OriginalData) {
foundOriginalData = true;
}
}
if (!foundOriginalData) {
checkBestBasis( top->lhsChild() );
}
if (!foundOriginalData) {
checkBestBasis( top->rhsChild() );
}
}
} // checkBestBasis
/**
Return true is be best basis has been calculated properly,
return false if the best basis has not been calculated or
it was not calculated properly.
The best basis is calculated in reference to a particular
cost function. A particular cost function will not always
result in a data set which differs from the original data
set. If this is the case, the "best basis" result will
include the original data.
*/
bool packtree::bestBasisOK()
{
foundOriginalData = false;
foundBestBasisVal = false;
checkBestBasis( root );
bool rslt = (foundBestBasisVal && (!foundOriginalData));
return rslt;
} // bestBasicOK
/**
Traverse the tree from the top down and add the best basis
nodes to the best basis list.
Note that the list object is simply a package for a scalar
value, the pointer to the head of the list. So it can be
passed by value without incurring a cost greater than
passing a pointer (e.g., pass by reference).
*/
void packtree::buildBestBasisList( packnode *top,
packdata_list &list )
{
if (top != 0) {
if (top->mark()) {
list.add( top );
}
else {
buildBestBasisList( top->lhsChild(), list );
buildBestBasisList( top->rhsChild(), list );
}
}
} // buildBestBasisList
/**
Return a list consisting of the best basis packdata values.
*/
packdata_list packtree::getBestBasisList()
{
packdata_list list;
buildBestBasisList( root, list );
return list;
} // getBestBasisList
packet/src/packtree_base.cpp 0100644 0000764 0000764 00000013131 07515163066 015057 0 ustar iank iank
#include "queue.h"
#include "packnode.h"
#include "packcontainer.h"
#include "packtree_base.h"
/**
Add a new level to the wavelet packet tree. Wavelet packet trees
are data structures that support a variety of applications.
If the reverse argument is true, the locations of the high pass and
low pass filter results in the wavelet calculation will be
reversed. This is used in building a wavelet packet tree for
frequency analysis.
The top packnode object contains data from the previous
level. If this is the first level in the tree, top will
contain the input data set.
A packcontainer object is created with the top data. The
packcontainer constructor will allocate new storage and copy the
data into this storage. If the packcontainer object consists of N
elements, then there will be N/2 elements on the left and N/2 or the
right. The object allows the lhs and rhs data to be accessed as one
array.
A wavelet transform step is calculated on the N element data set.
This results in N/2 values from the wavelet scaling function (low
pass function) and N/2 values from the wavelet function (high pass
function). These values are used to create two new packnode objects
which become children of top. Sub-trees for the new packnode
objects are recursively calculated.
As the tree is constructed, the leaves of the tree are marked with a
boolean flag in preparation for calculating the "best basis"
representation for the data. See the algorithm outlined in section
8.2.2 of "Ripples in Mathematics" by Jensen and la Cour-Harbo
The wavelet packet tree form that is used for wavelet packet
frequency analysis is described in section 9.3 (figure 9.14)
ofg "Ripples in Mathematics".
*/
void packtree_base::newLevel( packnode* top,
bool freqCalc,
bool reverse )
{
if (top != 0) {
const size_t len = top->length();
if (len > 1) {
// Create a new wavelet packet container for use in
// calculating the wavelet transform. Note that the
// container is only used locally.
packcontainer container( top );
if (reverse) {
// Calculate the reverse foward wavelet transform step,
// where the high pass result is stored in the upper half
// of the container and the low pass result is stored
// in the lower half of the container.
waveObj->forwardStepRev( container, len );
}
else {
// Calculate the foward wavelet transform step, where
// the high pass result is stored in the upper half
// of the container and the low pass result is stored
// in the lower half of the container.
waveObj->forwardStep( container, len );
}
packnode *lhs = new packnode(container.lhsData(),
len/2,
packnode::LowPass );
packnode *rhs = new packnode(container.rhsData(),
len/2,
packnode::HighPass );
// set the "mark" in the top node to false and
// mark the two children to true.
top->mark( false );
lhs->mark( true );
rhs->mark( true );
top->lhsChild( lhs );
top->rhsChild( rhs );
// The transform on the left hand side always uses
// the standard order (e.g., low pass filter result
// goes in the lower half, high pass goes in the
// upper half of the container).
newLevel( lhs, freqCalc, false );
if (freqCalc) {
// wavelet packet frequency analysis reverses the
// storage locations for the filter results in the
// right hand child
newLevel( rhs, freqCalc, true );
}
else { // freq == false
// use standard filter location
newLevel( rhs, freqCalc, false );
}
}
}
} // newLevel
/**
Print the wavelet packet tree, breadth first (this is also
sometimes called a level traversal).
The breadth first traversal uses a queue. The root of the tree is
initially inserted into the queue. The function operates by
deleting the node at teh front of the queue, printing the data
associated with that node and adding the node's left and right
children to the queue. Since a node's children are at the next
lower level, and we add the left child before the right child,
the function prints a tree level from left to right.
The above paragraph paraphrases Chapter 5, Level Order Traversal
of Fundamentals of Data Structures in C by Horowitz, Sahni
and Anderson-Freed, 1993.
*/
void packtree_base::breadthFirstPrint(const printKind kind)
{
queue Q;
Q.addQueue( root, 0 );
while (! Q.queueEmpty() ) {
packnode *node = Q.queueStart()->node;
size_t indent = Q.queueStart()->indent;
Q.deleteStart();
if (indent > 0) {
// print 'indent' spaces
printf("%*c", indent, ' ');
}
switch (kind) {
case printData: node->pr();
break;
case printCost: node->prCost();
break;
case printBestBasis: node->prBestBasis();
break;
default:
assert( false );
break;
} // switch
packnode *lhs = node->lhsChild();
packnode *rhs = node->rhsChild();
if (lhs != 0) {
Q.addQueue( lhs, indent + 2 );
}
if (rhs != 0) {
Q.addQueue( rhs, indent + 2 );
}
}
} // packtree_base::breadthFirstPrint
/**
Print the wavelet packet tree data and wavelet transform
result to standard out.
*/
void packtree_base::pr()
{
if (root != 0) {
breadthFirstPrint(printData);
}
} // pr
packet/src/blockpool.cpp~ 0100644 0000764 0000764 00000016156 07447470002 014456 0 ustar iank iank
/** \file
This file contains the class functions for a "pool" based memory
allocator. These functions manage the block chain that makes up the
memory pool and allocates memory from the individual blocks. The
entire pool is freed at once.
The pool allocator was originally written as a low level allocation
package. The original package queries the system to find the page
size and allocation block size. This version has been simplified,
but some of the low level nature remains.
The documentation in this file is formatted for doxygen
(see www.doxyeng.org).
Copyright and Use
You may use this source code without limitation and without
fee as long as you include:
This software was written and is copyrighted by Ian Kaplan, Bear
Products International, www.bearcave.com, 2002.
This software is provided "as is", without any warranty or
claim as to its usefulness. Anyone who uses this source code
uses it at their own risk. Nor is any support provided by
Ian Kaplan and Bear Products International.
Please send any bug fixes or suggested source changes to:
iank@bearcave.com
@author Ian Kaplan
*/
#include
#include
#include
#include "blockpool.h"
unsigned int block_pool::alloc_gran = (unsigned int)block_pool::page_size;
block_pool::block_chain *block_pool::current_block = 0;
block_pool::block_chain *block_pool::block_list_start = 0;
/**
block_pool::init_pool
This function is automatically called to initialize the block_pool
object when the "current_block" pointer is null.
*/
void block_pool::init_pool( void )
{
block_chain *new_link;
new_link = new_block( alloc_gran );
block_list_start = new_link;
current_block = new_link;
} /* init_pool */
/**
new_block
The new_block function is the "root" memory allocator for the
block_pool object. The amount of memory allocated is rounded up to
the next "block_size" boundary. Both the block_chain structure and
the allocatible memory are allocated from a single block that is a
multiple of the page size. This should avoid fragmentation in the
system memory allocator.
The "page" referenced here is a virtual memory page. This version
of the code does not actually check the system page size, but assumes
a 4Kb page.
*/
block_pool::block_chain *block_pool::new_block( unsigned int block_size )
{
const unsigned int max_block_size = max_block_multiple * page_size;
block_chain *new_link = 0;
unsigned int alloc_amt, total_alloc;
// add in the memory needed for the block_chain structure
total_alloc = block_size + sizeof(block_chain);
if (total_alloc < alloc_gran)
alloc_amt = alloc_gran;
else {
// its larger than the minimum allocation granularity, so round
// up the the nearest page.
alloc_amt = ((total_alloc + (page_size-1))/page_size) * page_size;
}
if (alloc_amt <= max_block_size) {
/* Allocate memory for both the block_chain structure and the memory
block */
new_link = (block_chain *)MemAlloc( alloc_amt );
// The new memory block starts after the block_chain structure
Chain_block(new_link) = (void *)(((unsigned int)new_link) + sizeof(block_chain));
assert( alloc_amt >= block_size );
Chain_bytes_used(new_link) = 0;
Chain_block_size(new_link) = alloc_amt - sizeof(block_chain);
Chain_next(new_link) = 0;
}
else {
printf("block_pool::new_block: allocation request too large\n");
}
return new_link;
} // block_chain::new_block
/**
block_pool::add_block Add a new memory block to the memory pool.
This function is called when the amount of memory requested by
pool_alloc will not fit in the current block.
*/
void *block_pool::add_block( unsigned int block_size )
{
block_chain *block = 0;
block_chain *last_block;
last_block = current_block;
block = new_block( block_size );
Chain_next(current_block) = block;
current_block = block;
return (void *)block;
} // block_chain::add_block
/**
pool_alloc
Allocate memory from the memory pool. The pool_alloc and free_pool
functions do memory allocation and deallocation.
This function is called to allocate memory from the memory pool.
If there is enough free memory in the current block to satisify the
memory request, memory is allocated from the current block and the
amount of free memory is updated. If the current block does not
have enough memory, add_block is called to allocate a new memory
block which will be large enough.
*/
void *block_pool::pool_alloc( unsigned int num_bytes )
{
const unsigned int align = sizeof( int );
void *addr = 0;
unsigned int amt_free;
/* the number of bytes allocated must be a multiple of the align
size */
num_bytes = ((num_bytes + (align-1))/align) * align;
if (current_block == 0) {
init_pool();
}
amt_free = Chain_block_size(current_block) - Chain_bytes_used(current_block);
if (num_bytes > amt_free) {
if (add_block( num_bytes ) != 0) {
amt_free = Chain_block_size(current_block);
}
}
if (amt_free >= num_bytes) {
addr = (void *)((unsigned int)Chain_block(current_block) + Chain_bytes_used(current_block));
Chain_bytes_used(current_block) += num_bytes;
}
else {
printf("block_pool::block_alloc: allocation error\n");
exit(1);
}
return addr;
} // block_pool::pool_alloc
/**
block_pool::free_pool
Walk through the block chain and deallocate the blocks.
Note that the block chain structures and the allocatible
memory is contained within a single allocated block.
The block_chain structure is at the start of this block
so passing its address to the memory deallocation function
deallocates both the block chain structure and the allocatible
memory
*/
void block_pool::free_pool(void)
{
block_chain *tmp;
while (block_list_start != 0) {
tmp = block_list_start;
block_list_start = Chain_next(block_list_start);
MemFree( (void *)tmp );
}
} // free_pool
/**
print_block_pool_info
Print information about the block pool
*/
void block_pool::print_block_pool_info( FILE *fp /*= stdout */)
{
int total_allocated = 0;
int total_unused = 0;
block_chain *ptr = block_list_start;
fprintf(fp, "Minimum memory allocation size: %d\n", alloc_gran );
fprintf(fp, "Page size: %d\n", (unsigned int)page_size );
fprintf(fp, "[block size, bytes_used]\n");
while (ptr != 0) {
fprintf(fp, "[%4d, %4d]", Chain_block_size(ptr), Chain_bytes_used(ptr));
total_allocated += Chain_bytes_used(ptr);
total_unused += (Chain_block_size(ptr) - Chain_bytes_used(ptr));
if (Chain_next(ptr) != 0) {
fprintf(fp, ", ");
}
else {
fprintf(fp, "\n");
}
ptr = Chain_next(ptr);
} // while
fprintf(fp, "Total allocated = %5d, total unused = %3d\n", total_allocated,
total_unused );
}