|
|
#include "complex.h"
Complex::Complex()
{
}
Complex::Complex(double r, double i)
: ComplexBase(r, i)
{}
void Complex::operator = (double scalar)
{
real(scalar);
imag(scalar);
}
const Complex Complex::operator / (double scalar)
{
Complex x(*this);
x /= scalar;
return x;
}
double Complex::length() const
{
return sqrt(real() * real() + imag() * imag());
}
double Complex::direction() const
{
return atan(imag() / real());
}
|
|
|
#pragma once
#include <complex>
using std::complex;
typedef complex<double> ComplexBase;
class Complex : public ComplexBase
{
public:
Complex();
Complex(double r, double i);
const Complex operator / (double scalar);
void operator = (double scalar);
double length() const;
double direction() const;
} ;
|
|
|
#include "complexvector.h"
void ComplexVector::operator /= (double scalar)
{
for (size_type i = 0; i < size(); i++)
at(i) /= scalar;
}
void ComplexVector::operator = (const ComplexVector& rhs)
{
for (size_type i = 0; i < size(); i++)
at(i) = rhs[i];
}
|
|
|
#pragma once
#include <vector>
using std::vector;
#include "Complex.h"
typedef vector<Complex> ComplexVectorBase;
class ComplexVector : public ComplexVectorBase
{
public:
void operator /= (double scalar);
void operator = (const ComplexVector& rhs);
} ;
|
|
|
#include "Signal.h"
#include "SignalGenerator.h"
int main()
{
Signal signal(128, 880*4);
SignalGenerator sgen(signal);
sgen.AddSine(7.0, 440.0, true);
sgen.AddSine(8.0, 880.0);
sgen.AddBias(2.0);
//signal.printTimeDomain("sampled data");
signal.ForwardDFT();
signal.printFreqDomain("result");
return 0;
}
|
|
|
#include <iostream>
#include <iomanip>
using namespace std;
#define _USE_MATH_DEFINES
#include <math.h>
#include "Signal.h"
Signal::Signal(long numsamples, long samplerate)
: mSampleRate(samplerate)
{
mData.resize(numsamples, Complex());
}
long Signal::NumSamples() //num samples
{
return (long) mData.size();
}
double Signal::SampleInterval()
{
return 1.0 / (double) mSampleRate;
}
void Signal::SetAll(ComplexFunc& f)
{
SetT(f);
}
void Signal::SetAll(DoubleFunc& f)
{
SetT(f);
}
void Signal::AddAll(ComplexFunc& f)
{
AddT(f);
}
void Signal::AddAll(DoubleFunc& f)
{
AddT(f);
}
static void fdbl( ios_base& os, double d)
{
dynamic_cast<ostream&>(os) << setw(8) << fixed << setprecision(3) << d;
}
static _Smanip<double> __cdecl fmtdouble(double d)
{
return (_Smanip<double>(&fdbl, d));
}
void Signal::printTimeDomain(const string& title)
{
cout << title << endl;
for(long i = 0; i < NumSamples(); ++i)
{
cout << setw(4) << (int) i
<< " "
<< fmtdouble(mData[i].real())
<< " "
<< fmtdouble(mData[i].imag())
<< endl;
}
}
void Signal::printFreqDomain(const string& title)
{
cout << title << endl;
cout << "Bias: ampl="
<< fmtdouble(mData[0].real())
<< endl;
double arg = (double) NumSamples() / (double)mSampleRate;
for(long i = 1; i < NumSamples() / 2 ; ++i)
{
double ampl = 2.0 * mData[i].length();
double phase = mData[i].direction();
cout << fmtdouble(i / arg)
<< " "
<< fmtdouble(ampl)
<< " "
<< fmtdouble(phase)
<< endl;
}
}
void Signal::ForwardDFT()
{
DFT(1);
}
void Signal::InverseDFT()
{
DFT(-1);
}
//direction = 1 => forward Transform
//direction = -1 => inverse Transform
void Signal::DFT(short int direction)
{
ComplexVector result;
result.resize(NumSamples());
double arg = -direction * 2.0 * M_PI / (double) NumSamples();
for (long f = 0; f < NumSamples(); f++)
for (long k = 0; k < NumSamples(); k++)
if (f == 0 || k == 0)
result[f] += mData[k] * Complex(1.0, 0.0);
else
result[f] += mData[k] * Complex(cos(k * f * arg), sin(k * f * arg));
if (direction == 1)
result /= (double) NumSamples();
mData = result;
}
|
|
|
#pragma once
#include <string>
using std::string;
#include "ComplexVector.h"
class Signal
{
public:
Signal(long numsamples, long samplerate);
void printTimeDomain(const string& title);
void printFreqDomain(const string& title);
void ForwardDFT();
void InverseDFT();
struct ComplexFunc
{
virtual Complex Value(int i, double t) = 0;
} ;
struct DoubleFunc
{
virtual double Value(int i, double t) = 0;
} ;
protected:
long NumSamples();
double SampleInterval();
void DFT(short int direction);
friend class SignalGenerator;
template <class T>
void AddT(T& f)
{
double t = 0.0;
for(long i = 0; i < NumSamples(); ++i, t += SampleInterval())
mData[i] += f.Value(i, t);
}
template <class T>
void SetT(T& f)
{
double t = 0.0;
for(long i = 0; i < NumSamples(); ++i, t += SampleInterval())
mData[i] = f.Value(i, t);
}
void SetAll(ComplexFunc& fn);
void AddAll(ComplexFunc& fn);
void SetAll(DoubleFunc& fn);
void AddAll(DoubleFunc& fn);
ComplexVector mData;
long mSampleRate;
} ;
|
|
|
#define _USE_MATH_DEFINES
#include <math.h>
#include "signalgenerator.h"
#include "complex.h"
#include "signal.h"
SignalGenerator::SignalGenerator(Signal& signal)
: mSignal(signal)
{
}
void SignalGenerator::AddSine(double ampl, double freq, bool clear)
{
struct Functor : public Signal::ComplexFunc
{
Functor(double a, double f) : ampl(a), freq(f) {}
Complex Value(int i, double t)
{
return Complex(t, ampl * sin(2.0 * M_PI * freq * t));
}
double ampl;
double freq;
} f(ampl, freq);
if (clear)
mSignal.SetAll(f);
else
mSignal.AddAll(f);
}
void SignalGenerator::AddBias(double ampl, bool clear)
{
struct Functor : public Signal::DoubleFunc
{
Functor(double a) : ampl(a) {}
double Value(int i, double t)
{
return ampl;
}
double ampl;
} f(ampl);
if (clear)
mSignal.SetAll(f);
else
mSignal.AddAll(f);
}
|
|
|
#pragma once
class Signal;
class SignalGenerator
{
public:
explicit SignalGenerator(Signal& signal);
void AddSine(double ampl, double freq, bool clear = false);
void AddBias(double ampl, bool clear = false);
private:
SignalGenerator();
SignalGenerator(const SignalGenerator& other);
Signal& mSignal;
} ;
|
|
|
#ifndef _BLOCKPOOL_H_
#define _BLOCKPOOL_H_
/** \file
The documentation in this file is formatted for doxygen
(see www.doxygen.org).
<b>Copyright and Use</b>
You may use this source code without limitation and without
fee as long as you include:
<blockquote>
This software was written and is copyrighted by Ian Kaplan, Bear
Products International, www.bearcave.com, 2002.
</blockquote>
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:
<pre>
iank@bearcave.com
</pre>
@author Ian Kaplan
*/
#include <assert.h>
#include <stdio.h>
#include <stdlib.h>
/**
This class supports memory pool allocation. A memory pool is
allocated in blocks and smaller chunks of memory are allocated from
these blocks. Instead of calling a set of class destructors (which
is time consuming) the memory pool allows all the objects allocated
in the memory pool to be destroyed at once. This provide a sort of
"poorman's garbage collector. Deallocating memory in a single place
simplifies the software structure, since allocation can be scattered
throughout the code without worry about deallocation.
This is a simplified version of a low level memory allocator that
can be used to create multiple memory pools. In this class
the class variables are static and shared by all instances of the
class. This allows the class to be declared locally to allocate
memory, but the state remains global. This makes the local
instance a window into the global state. The limitation is that
only one memory pool can be used.
Originally written November, 1996<br>
Revised for the wavelet packet transform code March 2002
\author Ian Kaplan
*/
class block_pool {
public: // typedefs and variables
/** the largest block of memory that can be allocated is the
page_size * max_block_multiple */
typedef enum { one_kay = 1024,
page_size = (4 * one_kay), /* 4 Kb */
max_block_multiple = 256, /* 1 Mb */
last_enum
} bogus;
/** typedef for memory block chain */
typedef struct block_chain_struct {
/** pointer to the current block */
void *block;
/** number of bytes used in the block */
unsigned int bytes_used;
/** total block size */
unsigned int block_size;
/** pointer to the next block */
block_chain_struct *next_block;
} block_chain;
private:
/** allocation granularity */
static unsigned int alloc_gran;
/** start of the block list for this pool */
static block_chain *block_list_start;
/** current block memory is being allocated from */
static block_chain *current_block;
private: // class functions
block_chain *new_block( unsigned int block_size );
void *add_block( unsigned int block_size );
void init_pool(void);
protected:
/**
Allocate memory using <i>calloc</i>. The POSIX calloc function
sets the memory to zero (in contrast to <i>malloc</i> which
allocates the memory without initializing it.
<pre>
#include <stdlib>
void *calloc( unsigned int num, unsigned int size );
</pre>
The arguments to calloc are:
<b>num</b>: number of elements<br>
<b>size</b>: size of the elements in bytes
The argument to MemAlloc is the size, in bytes
to allocate.
*/
virtual void *MemAlloc( unsigned int n_bytes )
{
void *rtn = calloc( n_bytes, 1 );
return rtn;
}
/**
Free memory that has been allocated with MemAlloc
*/
virtual void MemFree( void *addr )
{
free( addr );
}
public: // class functions
/** constructor does nothing, since all the class variables
are static */
block_pool(void) {}
void free_pool(void);
void *pool_alloc( unsigned int block_size );
void print_block_pool_info( FILE *fp = stdout );
}; // class block_pool
/* Macros for block_chain pointers */
/** return the memory block */
#define Chain_block(p) ((p)->block)
/** return the number of bytes used */
#define Chain_bytes_used(p) ((p)->bytes_used)
/** return the size of the memory block */
#define Chain_block_size(p) ((p)->block_size)
/** return the next block_chain structure in the list */
#define Chain_next(p) ((p)->next_block)
#endif
|
|
|
#ifndef _COSTBASE_H_
#define _COSTBASE_H_
#include "packnode.h"
/** \file
The documentation in this file is formatted for doxygen
(see www.doxygen.org).
<h4>
Copyright and Use
</h4>
You may use this source code without limitation and without
fee as long as you include:
<blockquote>
This software was written and is copyrighted by Ian Kaplan, Bear
Products International, www.bearcave.com, 2002.
</blockquote>
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:
<pre>
iank@bearcave.com
</pre>
@author Ian Kaplan
*/
/**
Base class for objects that define costs functions
for the wavelet packet transform.
The costbase base class provides a constructor that is passed the
root of a wavelet packet tree (which is built by packtree.cpp).
The 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 <i>Ripples in
Mathematics</i> by Jense and la Cour-Harbo.
\author Ian Kaplan
*/
class costbase
{
private:
/** disallow the copy constructor */
costbase( const costbase &rhs ) {}
protected:
/**
Recursively traverse the wavelet packet tree and calculate the cost
function.
*/
void traverse( packnode<double> *node )
{
if (node != 0) {
double cost = costCalc( node );
node->cost( cost );
traverse( node->lhsChild() );
traverse( node->rhsChild() );
}
} // traverse
/** Cost function to be defined by the subclass */
virtual double costCalc(packnode<double> *node) = 0;
public:
/** The default constructor does nothing */
costbase() {}
}; // costbase
#endif
|
|
|
#ifndef _COSTSHANNON_H_
#define _COSTSHANNON_H_
#include "costbase.h"
#include "packnode.h"
/** \file
The documentation in this file is formatted for doxygen
(see www.doxygen.org).
<b>Copyright and Use</b>
You may use this source code without limitation and without
fee as long as you include:
<blockquote>
This software was written and is copyrighted by Ian Kaplan, Bear
Products International, www.bearcave.com, 2002.
</blockquote>
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:
<pre>
iank@bearcave.com
</pre>
@author Ian Kaplan
*/
/**
The costshannon class extends the abstract class
costbase with a concrete implementation of the
costCalc function that implements the a modified
version of the Shannon entropy function as a cost
function.
\author Ian Kaplan
*/
class costshannon : public costbase
{
protected:
double costCalc( packnode<double> *node);
public:
/** Calculate a modified version of the the Shannon entropy cost
function for the wavelet packet tree, filling in the cost value
at each node */
costshannon( packnode<double> *node ) { traverse( node ); }
};
#endif
|
|
|
#ifndef _COSTTHRESH_H_
#define _COSTTHRESH_H_
#include "costbase.h"
#include "packnode.h"
/** \file
The documentation in this file is formatted for doxygen
(see www.doxygen.org).
<h4>
Copyright and Use
</h4>
<p>
You may use this source code without limitation and without
fee as long as you include:
</p>
<blockquote>
This software was written and is copyrighted by Ian Kaplan, Bear
Products International, www.bearcave.com, 2002.
</blockquote>
<p>
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.
<p>
Please send any bug fixes or suggested source changes to:
<pre>
iank@bearcave.com
</pre>
@author Ian Kaplan
*/
/**
The costthresh class implements a wavelet packet transform
cost function which counts the values in the data set that
are greater than a threshold value. This threshold value
is initialized in the constructor.
The costbase class provides the wavelet packet tree traversal
function which calls the subclass's cost function.
\author Ian Kaplan
*/
class costthresh : public costbase
{
private:
/** cost threshold */
double thresh;
/** absolute value */
double abs(const double x)
{
double retval = x;
if (retval < 0.0)
retval = -retval;
return retval;
} // abs
protected:
/**
This is a simple threshold calculation. The costCalc
function returns the number of node data values whose
absolute value is greater than the threshold.
*/
double costCalc(packnode<double> *node)
{
double count = 0.0;
if (node != 0) {
unsigned int len = node->length();
for (int i = 0; i < len; i++) {
if (abs((*node)[i]) > thresh) {
count = count + 1.0;
}
}
}
return count;
} // costCalc
public:
/** class constructor: calculate the wavelet packet cost
function using a simple threshold, t. */
costthresh(packnode<double> *node, double t )
{
thresh = t;
traverse( node );
}
}; // costthresh
#endif
|
|
|
#include <math.h>
/**
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 <i>Ripples in Mathematics</i> 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
<pre>
coef<sub>i</sub> = odd<sub>i</sub> - even<sub>i</sub>
or
coef<sub>i</sub> = 0.5 * (odd<sub>i</sub> - even<sub>i</sub>)
</pre>
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 h<sub>i</sub> and the
wavelet coefficients are g<sub>i</sub>.
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.
<pre>
a<sub>i</sub> = ...h0,h1,h2,h3, 0, 0, 0, 0, 0, 0, 0, ... s<sub>i</sub>
c<sub>i</sub> = ...g0,g1,g2,g3, 0, 0, 0, 0, 0, 0, 0, ... s<sub>i+1</sub>
a<sub>i+1</sub> = ...0, 0, h0,h1,h2,h3, 0, 0, 0, 0, 0, ... s<sub>i+2</sub>
c<sub>i+1</sub> = ...0, 0, g0,g1,g2,g3, 0, 0, 0, 0, 0, ... s<sub>i+3</sub>
a<sub>i+2</sub> = ...0, 0, 0, 0, h0,h1,h2,h3, 0, 0, 0, ... s<sub>i+4</sub>
c<sub>i+2</sub> = ...0, 0, 0, 0, g0,g1,g2,g3, 0, 0, 0, ... s<sub>i+5</sub>
a<sub>i+3</sub> = ...0, 0, 0, 0, 0, 0, h0,h1,h2,h3, 0, ... s<sub>i+6</sub>
c<sub>i+3</sub> = ...0, 0, 0, 0, 0, 0, g0,g1,g2,g3, 0, ... s<sub>i+7</sub>
</pre>
The dot product (inner product) of the infinite vector and
a row of the matrix produces either a smoother version of the
signal (a<sub>i</sub>) or a wavelet coefficient (c<sub>i</sub>).
In an ordered wavelet transform, the smoothed (a<sub>i</sub>) are
stored in the first half of an <i>n</i> element array region. The
wavelet coefficients (c<sub>i</sub>) are stored in the second half
the <i>n</i> 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.
<pre>
s<sub>i</sub> = ...h2,g2,h0,g0, 0, 0, 0, 0, 0, 0, 0, ... a<sub>i</sub>
s<sub>i+1</sub> = ...h3,g3,h1,g1, 0, 0, 0, 0, 0, 0, 0, ... c<sub>i</sub>
s<sub>i+2</sub> = ...0, 0, h2,g2,h0,g0, 0, 0, 0, 0, 0, ... a<sub>i+1</sub>
s<sub>i+3</sub> = ...0, 0, h3,g3,h1,g1, 0, 0, 0, 0, 0, ... c<sub>i+1</sub>
s<sub>i+4</sub> = ...0, 0, 0, 0, h2,g2,h0,g0, 0, 0, 0, ... a<sub>i+2</sub>
s<sub>i+5</sub> = ...0, 0, 0, 0, h3,g3,h1,g1, 0, 0, 0, ... c<sub>i+2</sub>
s<sub>i+6</sub> = ...0, 0, 0, 0, 0, 0, h2,g2,h0,g0, 0, ... a<sub>i+3</sub>
s<sub>i+7</sub> = ...0, 0, 0, 0, 0, 0, h3,g3,h1,g1, 0, ... c<sub>i+3</sub>
</pre>
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:
<pre>
a<sub>i</sub> = s[i]*h0 + s[i+1]*h1 + s[i+2]*h2 + s[i+3]*h3
c<sub>i</sub> = s[i]*g0 + s[i+1]*g1 + s[i+2]*g2 + s[i+3]*g3
</pre>
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 <i>Ripples in
Mathematics</i>, by Jensen and Cour-Harbo. Any error are mine.
<b>Author</b>: Ian Kaplan<br>
<b>Use</b>: 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 T>
class Daubechies : public liftbase<T, double> {
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
|
|
|
#ifndef FIFO_LIST_H
#define FIFO_LIST_H
/** \file
The documentation in this file is formatted for doxygen
(see www.doxygen.org).
<h4>
Copyright and Use
</h4>
<p>
You may use this source code without limitation and without
fee as long as you include:
</p>
<blockquote>
This software was written and is copyrighted by Ian Kaplan, Bear
Products International, www.bearcave.com, 2002.
</blockquote>
<p>
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.
<p>
Please send any bug fixes or suggested source changes to:
<pre>
iank@bearcave.com
</pre>
@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
<pre>
FIFO_LIST<char *> list;
FIFO_LIST<my_struct *> list;
FIFO_LIST<my_class *> list;
</pre>
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 T>
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(unsigned int 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
|
|
|
#ifndef _GROW_ARRAY_H_
#define _GROW_ARRAY_H_
/** \file
The documentation in this file is formatted for doxygen
(see www.doxyeng.org).
<h4>
Copyright and Use
</h4>
<p>
You may use this source code without limitation and without
fee as long as you include:
</p>
<blockquote>
This software was written and is copyrighted by Ian Kaplan, Bear
Products International, www.bearcave.com, 2002.
</blockquote>
<p>
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.
<p>
Please send any bug fixes or suggested source changes to:
<pre>
iank@bearcave.com
</pre>
@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 <vector> 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 T>
class GrowableArray {
private:
typedef enum { StartArraySize = 128 } bogus;
/** number of data elements */
unsigned int num_elem;
/** Array size (always <= num_elem) */
unsigned int 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;
unsigned int new_size = array_size * 2;
pArray = new T [ new_size ];
if (pArray != 0) {
rslt = true;
for (int i = 0; i < (int) 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 unsigned int length(void) const { return num_elem; }
/** set array to zero length */
void set_to_zero()
{
num_elem = 0;
}
/** LHS [] operator */
T &operator[](const unsigned int i)
{
assert( i < num_elem );
return pArray[ i ];
}
/** RHS [] operator */
T operator[](const unsigned int 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( unsigned int 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
|
|
|
#ifndef _HAAR_H_
#define _HAAR_H_
#include <math.h>
#include "liftbase.h"
/** \file
The documentation in this file is formatted for doxygen
(see www.doxygen.org).
<h4>
Copyright and Use
</h4>
<p>
You may use this source code without limitation and without
fee as long as you include:
</p>
<blockquote>
This software was written and is copyrighted by Ian Kaplan, Bear
Products International, www.bearcave.com, 2002.
</blockquote>
<p>
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.
<p>
Please send any bug fixes or suggested source changes to:
<pre>
iank@bearcave.com
</pre>
@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 even<sub>0</sub>, odd<sub>0</sub>,
even<sub>1</sub>, odd<sub>1</sub>, ...)
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 T>
class haar : public liftbase<T, double> {
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
<pre>
d<sub>j+1, i</sub> = odd<sub>j+1, i</sub> = odd<sub>j, i</sub> - even<sub>j, i</sub>
a<sub>j+1, i</sub> = even<sub>j, i</sub> = (even<sub>j, i</sub> + odd<sub>j, i</sub>)/2
</pre>
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
<pre>
a<sub>j+1, i</sub> = even<sub>j, i</sub> = even<sub>j, i</sub> + (odd<sub>j, i</sub>/2)
</pre>
*/
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
<pre>
double energy = 0.0;
for (int n = 0; n < N; n++) {
energy = energy + (a[i] * a[i]);
}
</pre>
See 5.2.1 of <i>Ripples in Mathematics</i> 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
|
|
|
#ifndef _HAAR_CLASSIC_H_
#define _HAAR_CLASSIC_H_
#include <math.h>
#include "liftbase.h"
/** \file
The documentation in this file is formatted for doxygen
(see www.doxygen.org)
<h4>
Copyright and Use
</h4>
<p>
You may use this source code without limitation and without
fee as long as you include:
</p>
<blockquote>
This software was written and is copyrighted by Ian Kaplan, Bear
Products International, www.bearcave.com, 2001.
</blockquote>
<p>
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.
<p>
Please send any bug fixes or suggested source changes to:
<pre>
iank@bearcave.com
</pre>
@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 <i>Ripples
in Mathematics</i> 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 a<sub>i</sub> is an even element
and an element b<sub>i</sub> is an odd element.
In this version of the Haar wavelet transform the scaling (or
smoothing) function is
<pre>
s = (a + b)/2
</pre>
The wavelet function is
<pre>
d = (a - b)/2
</pre>
A lifting scheme expression is used in this implementation. Here
the wavelet function is calculated first. The wavelet results
overwrite the odd b<sub>i</sub> values. This means that the
smoothing function values must be calculated with the result
of the wavelet function. To recover the value of b<sub>i</sub>
we use the expression
<pre>
b = a - 2d
s = (a + (a - 2d))/2
s = (2a - 2d)/2
s = a - d
</pre>
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 T>
class haar_classic : public liftbase<T, double> {
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
|
|
|
#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
<pre>
b'<sub>i</sub> = (a<sub>i</sub> - b<sub>i</i>)/2
</pre>
where a<sub>i</sub> is an even element and b<sub>i</sub> 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
<pre>
a'<sub>i</sub> = (a<sub>i</sub> - b<sub>i</i>)/2
</pre>
The "Haar classic" low pass filter (smoothing function) is
<pre>
a'<sub>i</sub> = (a<sub>i</sub> + b<sub>i</sub>)/2
</pre>
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
<pre>
b'<sub>i</sub> = (a<sub>i</sub> + b<sub>i</sub>)/2
</pre>
The calculation of the high pass filter is done first and
overwrites the even elements (e.g., a'<sub>i</sub>). To
recover the a<sub>i</sub> values given a'<sub>i</sub> and
b<sub>i</sub>:
<pre>
a'<sub>i</sub> = (a<sub>i</sub> - b<sub>i</i>)/2
2 * a'<sub>i</sub> = a<sub>i</sub> - b<sub>i</i>
a<sub>i</sub> = 2 * a'<sub>i</sub> + b<sub>i</i>
</pre>
Substituting this into the equation for calculating
b'<sub>i</i> we get:
<pre>
b'<sub>i</sub> = ((2 * a'<sub>i</sub> + b<sub>i</i>) + b<sub>i</sub>)/2
b'<sub>i</sub> = (2 * a'<sub>i</sub> + 2 * b<sub>i</sub>)/2
b'<sub>i</sub> = a'<sub>i</sub> + b<sub>i</sub>
</pre>
These equations differ from the Haar classic forward transform
equations.
\author Ian Kaplan
*/
template <class T>
class haar_classicFreq : public haar_classic<T> {
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
|
|
|
#ifndef _INVPACKTREE_H_
#define _INVPACKTREE_H_
/** \file
The documentation in this file is formatted for doxygen
(see www.doxygen.org).
<h4>
Copyright and Use
</h4>
<p>
You may use this source code without limitation and without
fee as long as you include:
</p>
<blockquote>
This software was written and is copyrighted by Ian Kaplan, Bear
Products International, www.bearcave.com, 2002.
</blockquote>
<p>
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.
<p>
Please send any bug fixes or suggested source changes to:
<pre>
iank@bearcave.com
</pre>
@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
<i>Ripples in Mathematics</i> 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 <i>Ripples in
Mathematics</i>.
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<packcontainer, double> *waveObj;
/** disallow the copy constructor */
invpacktree( const invpacktree &rhs ) {}
/** inverse wavelet packet transform calculation stack */
LIST<packcontainer *> stack;
/** pointer to the inverse transform result */
const double *data;
/** length of data */
unsigned int N;
private:
void new_level( packdata<double> *elem );
void add_elem( packdata<double> *elem );
void reduce();
public:
invpacktree( packdata_list<double> &list,
liftbase<packcontainer, double> *w );
/** The destructor does nothing */
~invpacktree() {}
/** Get the result of the inverse packet transform */
const double *getData() { return data; }
void pr();
};
#endif
|
|
|
#ifndef _LIFTBASE_H_
#define _LIFTBASE_H_
/** \file
<b>Copyright and Use</b>
You may use this source code without limitation and without
fee as long as you include:
<blockquote>
This software was written and is copyrighted by Ian Kaplan, Bear
Products International, www.bearcave.com, 2002.
</blockquote>
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:
<pre>
iank@bearcave.com
</pre>
@author Ian Kaplan
*/
#include <assert.h>
/**
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:
<ul>
<li>
<p>
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.
</p>
</li>
<li>
<p>
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.
</p>
</li>
<li>
<p>
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.
</p>
</li>
<li>
<p>
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.
</p>
</li>
</ul>
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.
<b>Instantiating the Template</b>
The liftbase template takes two type arguments:
<ol>
<li>
The type of the array or '[]' operator indexable object.
</li>
<li>
The type of the data element.
</li>
</ol>
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:
<pre>
class Haar : public liftbase<double *, double>
</pre>
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.
<pre>
class Poly : public liftbase<packcontainer, double>
</pre>
<b>References:</b>
<ul>
<li>
<a href="http://www.bearcave.com/misl/misl_tech/wavelets/packet/index.html">
<i>The Wavelet Packet Transform</i></a> by Ian Kaplan, www.bearcave.com.
</li>
<li>
<a
href="http://www.bearcave.com/misl/misl_tech/wavelets/lifting/index.html">
<i>The Wavelet Lifting Scheme</i></a> by Ian Kaplan, www.bearcave.com.
This is the parent web page for this Java source code.
</li>
<li>
<i>Ripples in Mathematics: the Discrete Wavelet Transform</i>
by Arne Jense and Anders la Cour-Harbo, Springer, 2001
</li>
<li>
<i>Building Your Own Wavelets at Home</i> in <a
href="http://www.multires.caltech.edu/teaching/courses/waveletcourse/">
Wavelets in Computer Graphics</a>
</li>
</ul>
\author Ian Kaplan
*/
template <class T, class T_elem >
class liftbase {
protected:
typedef enum {
/** "enumeration" for forward wavelet transform */
forward = 1,
/** "enumeration" for inverse wavelet transform */
inverse = 2
} transDirection;
/**
Split the <i>vec</i> 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
|
|
|
#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.
<h4>
Copyright and Use
</h4>
You may use this source code without limitation and without
fee as long as you include:
<blockquote>
This software was written and is copyrighted by Ian Kaplan, Bear
Products International, www.bearcave.com, 2001.
</blockquote>
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:
<pre>
iank@bearcave.com
</pre>
\author Ian Kaplan
*/
template <class T>
class line : public liftbase<T, double> {
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, x<sub>1</sub>, y<sub>1</sub> and
x<sub>2</sub>, y<sub>2</sub>, where
<pre>
x<sub>1</sub> = 0
x<sub>2</sub> = 1
</pre>
calculate the point on the line at x<sub>3</sub>, y<sub>3</sub>,
where
<pre>
x<sub>3</sub> = 2
</pre>
The "two-point equation" for a line given x<sub>1</sub>,
y<sub>1</sub> and x<sub>2</sub>, y<sub>2</sub> is
<pre>
. y<sub>2</sub> - y<sub>1</sub>
(y - y<sub>1</sub>) = -------- (x - x<sub>1</sub>)
. x<sub>2</sub> - x<sub>1</sub>
</pre>
Solving for y
<pre>
. y<sub>2</sub> - y<sub>1</sub>
y = -------- (x - x<sub>1</sub>) + y<sub>1</sub>
. x<sub>2</sub> - x<sub>1</sub>
</pre>
Since x<sub>1</sub> = 0 and x<sub>2</sub> = 1
<pre>
. y<sub>2</sub> - y<sub>1</sub>
y = -------- (x - 0) + y<sub>1</sub>
. 1 - 0
</pre>
or
<pre>
y = (y<sub>2</sub> - y<sub>1</sub>)*x + y<sub>1</sub>
</pre>
We're calculating the value at x<sub>3</sub> = 2, so
<pre>
y = 2*y<sub>2</sub> - 2*y<sub>1</sub> + y<sub>1</sub>
</pre>
or
<pre>
y = 2*y<sub>2</sub> - y<sub>1</sub>
</pre>
*/
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 (b<sub>j</sub>) to the second half
of the array, leaving the even elements (a<sub>i</sub>) in the
first half
<pre>
a<sub>0</sub>, a<sub>1</sub>, a<sub>1</sub>, a<sub>3</sub>, b<sub>0</sub>, b<sub>1</sub>, b<sub>2</sub>, b<sub>2</sub>,
</pre>
The predict step of the line wavelet "predicts" that the
odd element will be on a line between two even elements.
<pre>
b<sub>j+1,i</sub> = b<sub>j,i</sub> - (a<sub>j,i</sub> + a<sub>j,i+1</sub>)/2
</pre>
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 2<sup>n</sup>
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 <i>new_y()</i>
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
(<i>Building Your Own Wavelets at Home</i>) and in
<i>Ripples in Mathematics</i>. A somewhat more complete derivation
of this equation is provided in <i>Ripples in Mathematics</i>
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>i</i>. Note that the predict phase has already
completed, so the odd values belong to iteration <i>i+1</i>.
<pre>
even<sub>i+1,j</sub> = even<sub>i,j</sub> op (odd<sub>i+1,k-1</sub> + odd<sub>i+1,k</sub>)/4
</pre>
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 odd<sub>i+1,k-1</sub> is
the same as odd<sub>k</sub>. So for the first element this
becomes
<pre>
(2 * odd<sub>k</sub>)/4
</pre>
or
<pre>
odd<sub>k</sub>/2
</pre>
*/
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
|
|
|
#ifndef _LIST_H_
#define _LIST_H_
/** \file
The documentation in this file is formatted for doxygen
(see www.doxygen.org).
<h4>
Copyright and Use
</h4>
<p>
You may use this source code without limitation and without
fee as long as you include:
</p>
<blockquote>
This software was written and is copyrighted by Ian Kaplan, Bear
Products International, www.bearcave.com, 2002.
</blockquote>
<p>
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.
<p>
Please send any bug fixes or suggested source changes to:
<pre>
iank@bearcave.com
</pre>
@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
<pre>
LIST<char *> list;
LIST<my_struct *> list;
</pre>
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 T>
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(unsigned int 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<T> &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
|
|
|
#ifndef _PACKCONTAINER_H_
#define _PACKCONTAINER_H_
#include "packnode.h"
/**
The documentation in this file is formatted for doxygen
(see www.doxygen.org).
<h4>
Copyright and Use
</h4>
<p>
You may use this source code without limitation and without
fee as long as you include:
</p>
<blockquote>
This software was written and is copyrighted by Ian Kaplan, Bear
Products International, www.bearcave.com, 2002.
</blockquote>
<p>
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.
<p>
Please send any bug fixes or suggested source changes to:
<pre>
iank@bearcave.com
</pre>
@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 */
unsigned int 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<double>* node )
{
assert( node != 0 );
N = node->length();
assert( N > 1 );
unsigned int half = N >> 1;
block_pool mem_pool;
unsigned int num_bytes = half * sizeof(double);
lhs = (double *)mem_pool.pool_alloc( num_bytes );
rhs = (double *)mem_pool.pool_alloc( num_bytes );
for (unsigned int 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( unsigned int 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 unsigned int i )
{
assert( i < N );
unsigned int half = N >> 1;
if (i < half)
return lhs[i];
else {
return rhs[i-half];
}
}
/** RHS [] operator */
double operator[]( const unsigned int i ) const
{
assert( i < N );
unsigned int 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.
*/
unsigned int length() { return N; }
}; // packcontainer
#endif
|
|
|
#ifndef _PACKDATA_H_
#define _PACKDATA_H_
/** \file
The documentation in this file is formatted for doxygen
(see www.doxygen.org).
<h4>
Copyright and Use
</h4>
<p>
You may use this source code without limitation and without
fee as long as you include:
</p>
<blockquote>
This software was written and is copyrighted by Ian Kaplan, Bear
Products International, www.bearcave.com, 2002.
</blockquote>
<p>
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.
<p>
Please send any bug fixes or suggested source changes to:
<pre>
iank@bearcave.com
</pre>
@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 T>
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 */
unsigned int N;
/** Wavelet packet data */
T* data;
/** default constructor */
packdata() {}
public:
/**
\arg vec a pointer to an array of type <i>T</i> elements
\arg n the size of the array of type <i>T</i> elements
\arg k the kind of data (e.g., original data, low pass result,
high pass result.
*/
packdata( T *vec, const unsigned int n, const transformKind k )
{
data = vec;
N = n;
kind = k;
}
/** define a virtual destructure which in the base class */
virtual ~packdata() {}
/**
Overload the standard <i>new</i> 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 < (int) 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 */
unsigned int length() { return N; }
}; // packdata
#endif
|
|
|
#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).
<h4>
Copyright and Use
</h4>
<p>
You may use this source code without limitation and without
fee as long as you include:
</p>
<blockquote>
This software was written and is copyrighted by Ian Kaplan, Bear
Products International, www.bearcave.com, 2002.
</blockquote>
<p>
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.
<p>
Please send any bug fixes or suggested source changes to:
<pre>
iank@bearcave.com
</pre>
\author Ian Kaplan
*/
/**
The packdata_list class is based on a subclass built from a
FIFO_LIST<packdata *> 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 T>
class packdata_list : public FIFO_LIST<packdata<T> *>
{
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<T> *elem;
for (h = first(); h != 0; h = next( h )) {
elem = get_item( h );
elem->pr();
} // for h
printf("\n");
} // pr
}; // packdata_list
#endif
|
|
|
#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).
<h4>
Copyright and Use
</h4>
<p>
You may use this source code without limitation and without
fee as long as you include:
</p>
<blockquote>
This software was written and is copyrighted by Ian Kaplan, Bear
Products International, www.bearcave.com, 2002.
</blockquote>
<p>
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.
<p>
Please send any bug fixes or suggested source changes to:
<pre>
iank@bearcave.com
</pre>
\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:
<ol>
<li>
A vector of doubles containing the data set (e.g., the signal).
The length of the vector must be a power of two.
</li>
<li>
A pointer to a wavelet lifting scheme class that will be
used in calculating the wavelet transform step.
</li>
</ol>
If the vector passed to the constructor contains N double values, the
result of the constructor will be a wavelet packet tree with
log<sub>2</sub>(N) levels.
\author Ian Kaplan
*/
class packfreq : public packtree_base {
private:
/** Level basis matrix */
GrowableArray<packnode<double> *> mat;
void findLevel( packnode<double>* top,
unsigned int cur_level,
const unsigned int level );
protected:
/** disallow the copy constructor */
packfreq( const packfreq &rhs ) {};
/** disallow the default constructor */
packfreq() {};
public:
packfreq( const double *vec,
const unsigned int n,
liftbase<packcontainer, double> *w );
/** destructor does nothing */
~packfreq() {}
void getLevel( const unsigned int level );
void plotMat(const unsigned int N);
void prMat();
}; // packfreq
#endif
|
|
|
#ifndef _PACKNODE_H_
#define _PACKNODE_H_
#include "packdata.h"
/** \file
The documentation in this file is formatted for doxygen
(see www.doxygen.org).
<h4>
Copyright and Use
</h4>
<p>
You may use this source code without limitation and without
fee as long as you include:
</p>
<blockquote>
This software was written and is copyrighted by Ian Kaplan, Bear
Products International, www.bearcave.com, 2002.
</blockquote>
<p>
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.
<p>
Please send any bug fixes or suggested source changes to:
<pre>
iank@bearcave.com
</pre>
@author Ian Kaplan
*/
/**
A wavelet packet tree node
A description of the wavelet packet algorithm can be found
in Chapter 8 of <i>Ripples in Mathematics</i> 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 log<sub>2</sub>(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 <i>n</i> 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 T>
class packnode : public packdata<T> {
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 unsigned int n,
const transformKind k ) : packdata<T>(vec, n, k)
{
leftChild = 0;
rightChild = 0;
costVal = 0.0;
chosen = false;
}
/** LHS [] operator */
T &operator[]( const unsigned int i )
{
assert( i < N );
return data[i];
}
/** RHS [] operator */
T operator[]( const unsigned int 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 < (int) 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
|
|
|
#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).
<h4>
Copyright and Use
</h4>
<p>
You may use this source code without limitation and without
fee as long as you include:
</p>
<blockquote>
This software was written and is copyrighted by Ian Kaplan, Bear
Products International, www.bearcave.com, 2002.
</blockquote>
<p>
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.
<p>
Please send any bug fixes or suggested source changes to:
<pre>
iank@bearcave.com
</pre>
@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 log<sub>2</sub>(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<double> *root );
void buildBestBasisList( packnode<double> *root,
packdata_list<double> &list );
void checkBestBasis( packnode<double> *root );
void cleanTree(packnode<double> *root, bool removeMark );
public:
packtree( const double *vec,
const unsigned int n,
liftbase<packcontainer, double> *w );
/** destructor does nothing */
~packtree() {}
void prCost();
void prBestBasis();
void bestBasis();
bool bestBasisOK();
packdata_list<double> getBestBasisList();
}; // packtree
#endif
|
|
|
#ifndef _PACKTREE_BASE_H_
#define _PACKTREE_BASE_H_
/** \file
The documentation in this file is formatted for doxygen
(see www.doxygen.org).
<h4>
Copyright and Use
</h4>
<p>
You may use this source code without limitation and without
fee as long as you include:
</p>
<blockquote>
This software was written and is copyrighted by Ian Kaplan, Bear
Products International, www.bearcave.com, 2002.
</blockquote>
<p>
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.
<p>
Please send any bug fixes or suggested source changes to:
<pre>
iank@bearcave.com
</pre>
\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<double> *root;
/** wavelet packet transform object */
liftbase<packcontainer, double> *waveObj;
typedef enum { BadPrintKind,
printData,
printCost,
printBestBasis } printKind;
void breadthFirstPrint(printKind kind);
void newLevel( packnode<double>* top, bool freqCalc, bool reverse );
public:
void pr();
/** get the root of the wavelet packet tree */
packnode<double> *getRoot() { return root; }
}; // packtree_base
#endif
|
|
|
#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 <i>Ripples in Mathematics</i> 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).
<h4>
Copyright and Use
</h4>
<p>
You may use this source code without limitation and without
fee as long as you include:
</p>
<blockquote>
This software was written and is copyrighted by Ian Kaplan, Bear
Products International, www.bearcave.com, 2002.
</blockquote>
<p>
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.
<p>
Please send any bug fixes or suggested source changes to:
<pre>
iank@bearcave.com
</pre>
@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 T>
class queueElem
{
private:
/** disallow the default constructor */
queueElem() {}
public:
packnode<T> *node;
unsigned int indent;
/** Initialize the queue element */
queueElem( packnode<T> *n, unsigned int i )
{
node = n;
indent = i;
}
/** allocate queueElem objects from a 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
}; // 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 T>
class queue : protected FIFO_LIST<queueElem<T> *>
{
public:
/** Get the first element in the queue */
queueElem<T> *queueStart()
{
handle h = first();
queueElem<T> *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<T> *elem = get_item( h );
remove();
// no delete elem;
}
} // deleteStart
/** Add an element to the queue */
void addQueue(packnode<T> *node, unsigned int indent )
{
queueElem<T> *elem = new queueElem<T>(node, indent);
add( elem );
} // addQueue
/** return true if the queue is empty */
bool queueEmpty() { return (first() == 0); }
}; // queue
#endif
|
|
|
/** \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).
<h4>
Copyright and Use
</h4>
<p>
You may use this source code without limitation and without
fee as long as you include:
</p>
<blockquote>
This software was written and is copyrighted by Ian Kaplan, Bear
Products International, www.bearcave.com, 2002.
</blockquote>
<p>
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.
<p>
Please send any bug fixes or suggested source changes to:
<pre>
iank@bearcave.com
</pre>
@author Ian Kaplan
*/
#include <assert.h>
#include <stdio.h>
#include <stdlib.h>
#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 long)(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 long)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 );
}
|
|
|
#include <assert.h>
#include <stdio.h>
#include <math.h>
#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 <i>Ripples in Mathematics</i> 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<double> *node )
{
assert( node != 0 );
unsigned int len = node->length();
const double *a = node->getData();
double sum = 0.0;
for (int i = 0; i < (int) 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
|
|
|
/** \file
This file contains code to test wavelet frequency analysis via
wavelet packets. See
<a href="http://www.bearcave.com/misl/misl_tech/wavelets/packfreq/index.html">http://www.bearcave.com/misl/misl_tech/wavelets/packfreq/index.html</a>
The documentation in this file is formatted for doxygen
(see www.doxygen.org).
<h4>
Copyright and Use
</h4>
<p>
You may use this source code without limitation and without
fee as long as you include:
</p>
<blockquote>
This software was written and is copyrighted by Ian Kaplan, Bear
Products International, www.bearcave.com, 2002.
</blockquote>
<p>
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.
<p>
Please send any bug fixes or suggested source changes to:
<pre>
iank@bearcave.com
</pre>
@author Ian Kaplan
*/
#include <assert.h>
#include <stdio.h>
#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, unsigned int N)
{
const double PI = 3.1415926535897932384626433832795;
const double omega = 2 * PI;
const double samplerate = 1000; //N / 5.0; //128.0; //samples per second
const double time_incr = 1.0 / samplerate;
double t = 0.0;
for (int i = 0; i < (int) N; i++) {
vecX[i] = t;
vecY[i] = 0.0;
double ampl = 2.0; double freq = 3.0;
vecY[i] += ampl * sin(freq * omega * t );
//vecY[i] += 6 * sin( 441.0 * point );
//vecY[i] += 3 * sin( 220.0 * point );
//vecY[i] += 1 * sin( 880.0 * point );
//vecY[i] += 4 * sin( 64 * point );
//vecY[i] += 2 * sin( 32 * point );
//vecY[i] += 1 * sin( 16 * point );
//vecY[i] += 0.5 * sin( 8 * point );
// printf("%7.4g, %7.4g\n", vecX[i], vecY[i] );
t += time_incr;
}
} // gen_freqMix
/** \function
Generate a signal composed of one or more sine waves.
*/
void gen_sinCombo( double *vecX, double *vecY, unsigned int 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 < (int) 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, unsigned int N, unsigned int 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 unsigned int stepWidth = N / steps;
int stepCnt = 0;
int i;
for (i = 0; i < (int) 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, unsigned int 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 < (int) 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, unsigned int len )
{
for (int i = 0; i < (int) len; i++) {
printf("%7.4f %7.4f\n", vecX[i], vecY[i] );
}
}
/** \function
Print a vector of doubles
*/
void prVec( double *vec, unsigned int len )
{
for (int i = 0; i < (int) 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 unsigned int N = 1024;
// const unsigned int 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<packcontainer> h;
// Daubechies<packcontainer> 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.getLevel( 10);
//tree.prMat();
tree.plotMat(N);
// free the memory pool
block_pool mem_pool;
mem_pool.free_pool();
return 0;
}
|
|
|
#include <assert.h>
#include <stdio.h>
#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).
<h4>
Copyright and Use
</h4>
<p>
You may use this source code without limitation and without
fee as long as you include:
</p>
<blockquote>
This software was written and is copyrighted by Ian Kaplan, Bear
Products International, www.bearcave.com, 2002.
</blockquote>
<p>
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.
<p>
Please send any bug fixes or suggested source changes to:
<pre>
iank@bearcave.com
</pre>
@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<double> *elem )
{
unsigned int half = elem->length();
unsigned int 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<packcontainer *>::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();
unsigned int 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 < (int) 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<double> *elem )
{
assert( elem != 0 );
if (stack.first() == 0) {
new_level( elem );
}
else {
unsigned int half = elem->length();
unsigned int n = half * 2;
LIST<packcontainer *>::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<double> &list,
liftbase<packcontainer, double> *w )
{
data = 0;
N = 0;
waveObj = w;
// Traverse the "best basis" list and calculate the inverse
// wavelet packet transform.
packdata_list<double>::handle h;
for (h = list.first(); h != 0; h = list.next( h )) {
packdata<double> *elem = list.get_item( h );
add_elem( elem );
} // for
LIST<packcontainer *>::handle tosHandle;
tosHandle = stack.first();
packcontainer *tos = stack.get_item( tosHandle );
if (tos != 0) {
unsigned int 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 < (int) N; i++) {
printf("%7.4f ", data[i] );
}
printf("\n");
}
} // pr
|
|
|
#include <stdio.h>
#include <stdlib.h>
/** \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( unsigned int num_bytes )
{
printf("global operator new\n");
void *rtn = malloc( num_bytes );
return rtn;
} // new
void *operator new[]( unsigned int 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 );
}
|
|
|
/** \file
The documentation in this file is formatted for doxygen
(see www.doxyeng.org).
<h4>
Copyright and Use
</h4>
<p>
You may use this source code without limitation and without
fee as long as you include:
</p>
<blockquote>
This software was written and is copyrighted by Ian Kaplan, Bear
Products International, www.bearcave.com, 2002.
</blockquote>
<p>
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.
<p>
Please send any bug fixes or suggested source changes to:
<pre>
iank@bearcave.com
</pre>
@author Ian Kaplan
*/
#include <assert.h>
#include <stdio.h>
#include <math.h>
#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
<i>w</i>. 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 <b>H</b> and <b>G</b> 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 <b>H</b> and <b>G</b> results are in the upper and
lower halves of the result array. However, the left child inverts
the locations of the <b>H</b> and <b>G</b> results, placing the
result of the <b>G</b> filter in the lower half and the result
of the <b>H</b> 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 unsigned int N,
liftbase<packcontainer, double> *w )
{
waveObj = w;
block_pool mem_pool;
double *vecCopy = (double *)mem_pool.pool_alloc( N * sizeof( double ) );
for (int i = 0; i < (int) N; i++) {
vecCopy[i] = vec[i];
}
root = new packnode<double>( vecCopy, N, packnode<double>::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 <i>level</i> to the level basis matrix.
Levels are numbered from 0 at the root.
*/
void packfreq::findLevel( packnode<double>* top,
unsigned int cur_level,
const unsigned int 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.
<pre>
<b>number of elements</b>
<b>in a tree node</b> <b>level</b> <b>number of nodes</b>
1024 0 1
512 1 2
256 2 4
128 3 8
64 4 16
32 5 32
</pre>
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 unsigned int 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 <i>Ripples in Mathematics</i> by Jensen and
la-Cour Harbo (Springer Verlag, 2001), the calculation for the
z-axis value is:
<pre>
m = (*mat)[y][x]
z_val = log( 1 + val^2 )
</pre>
Where <i>log</i> 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 unsigned int N)
{
unsigned int num_y = mat.length();
const double incr = (double)N / (double)num_y;
if (num_y > 0) {
unsigned int num_x = mat[0]->length();
double freq_start = 0.0;
for (unsigned int y = 0; y < num_y; y++) {
double time_start = 0.0;
for (unsigned int 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(" %3.3d %3.3d %7.4f\n", y, x, log(1+(val*val)) );
// plot actual frequency and time values,
// with frequency on y and time on x
printf(" %9.4f %9.4f %9.4f\n",
time_start, freq_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()
{
unsigned int num_y = mat.length();
if (num_y > 0) {
unsigned int num_x = mat[0]->length();
for (int y = num_y-1; y >= 0; y--) {
for (unsigned int x = 0; x < num_x; x++) {
printf(" %7.4f ", (*mat[y])[ x ] );
}
printf("\n");
fflush(stdout);
}
}
} // prMat
|
|
|
/** \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).
<h4>
Copyright and Use
</h4>
<p>
You may use this source code without limitation and without
fee as long as you include:
</p>
<blockquote>
This software was written and is copyrighted by Ian Kaplan, Bear
Products International, www.bearcave.com, 2002.
</blockquote>
<p>
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.
<p>
Please send any bug fixes or suggested source changes to:
<pre>
iank@bearcave.com
</pre>
@author Ian Kaplan
*/
#include <assert.h>
#include <stdio.h>
#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<double *> 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<packcontainer> 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<double> *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<double> 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;
}
|
|
|
/** \file
The documentation in this file is formatted for doxygen
(see www.doxyeng.org).
<h4>
Copyright and Use
</h4>
<p>
You may use this source code without limitation and without
fee as long as you include:
</p>
<blockquote>
This software was written and is copyrighted by Ian Kaplan, Bear
Products International, www.bearcave.com, 2002.
</blockquote>
<p>
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.
<p>
Please send any bug fixes or suggested source changes to:
<pre>
iank@bearcave.com
</pre>
@author Ian Kaplan
*/
#include <assert.h>
#include <stdio.h>
#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 <i>w</i>. 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 unsigned int N,
liftbase<packcontainer, double> *w )
{
waveObj = w;
block_pool mem_pool;
double *vecCopy = (double *)mem_pool.pool_alloc( N * sizeof( double ) );
for (int i = 0; i < (int) N; i++) {
vecCopy[i] = vec[i];
}
root = new packnode<double>( vecCopy, N, packnode<double>::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 <i>Ripples in Mathematics</i>
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<double> *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 <i>Ripples in Mathematics</i>. 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<double> *top )
{
double cost = 0.0;
if (top != 0) {
packnode<double> *lhs = top->lhsChild();
packnode<double> *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<double> *top )
{
if (top != 0) {
if (top->mark()) {
foundBestBasisVal = true;
if (top->getKind() == packdata<double>::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<double> *top,
packdata_list<double> &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<double> packtree::getBestBasisList()
{
packdata_list<double> list;
buildBestBasisList( root, list );
return list;
} // getBestBasisList
|
|
|
#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 <i>top</i> packnode object contains data from the previous
level. If this is the first level in the tree, <i>top</i> will
contain the input data set.
A packcontainer object is created with the <i>top</i> 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 <i>top</i>. 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<double>* top,
bool freqCalc,
bool reverse )
{
if (top != 0) {
const unsigned int 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<double> *lhs = new packnode<double>(container.lhsData(),
len/2,
packnode<double>::LowPass );
packnode<double> *rhs = new packnode<double>(container.rhsData(),
len/2,
packnode<double>::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 <i>Fundamentals of Data Structures in C</i> by Horowitz, Sahni
and Anderson-Freed, 1993.
*/
void packtree_base::breadthFirstPrint(const printKind kind)
{
queue<double> Q;
Q.addQueue( root, 0 );
while (! Q.queueEmpty() ) {
packnode<double> *node = Q.queueStart()->node;
unsigned int 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<double> *lhs = node->lhsChild();
packnode<double> *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
|
|
|
#ifndef _YAHOOTS_H_
#define _YAHOOTS_H_
/** \file
The documentation in this file is formatted for doxygen
(see www.doxygen.org).
<h4>
Copyright and Use
</h4>
<p>
You may use this source code without limitation and without
fee as long as you include:
</p>
<blockquote>
This software was written and is copyrighted by Ian Kaplan, Bear
Products International, www.bearcave.com, 2002.
</blockquote>
<p>
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.
<p>
Please send any bug fixes or suggested source changes to:
<pre>
iank@bearcave.com
</pre>
@author Ian Kaplan
*/
/** Process historical equity (stock) data downloaded from
finance.yahoo.com. The data is downloaded in "spread sheet" format
from the historical data page. There is probably some limitation on
using this data (e.g., no commercial use and no resale) so use at
your own risk.
The format of the file is ASCII. The first line lists the title
for each of the fields in the file. The titles and the fields
are comma separated.
This class is specific to the data format that was downloaded
from Yahoo at the time. More general code could be written to
easily account for changing formats. However, I just wanted
to extract the data.
The Yahoo data has two places of accuracy, presumably reflecting
decimalization. The equity time series are adjusted for splits
and dividends, backward in time from the most recent time in the
time series. This can cause problems over long periods of time
since at some point a stock that pays dividends will pay all of
it's worth out in dividends and the value will become negative
(as a result, a reinvest is a better choice).
The format for the data is:
<title line>
<time series line>+
(e.g,. a titled followed by one or more time series lines).
The title line consists of six comma separated strings
(e.g., "Date,Open,High,Low,Close,Volume"). Time time
series lines have the values suggested in the title.
For my current purposes I am not interested in date values,
so these are ignored. All values are returned as vectors
of doubles, although volume is an unsigned integer value.
*/
class yahooTS
{
private:
const char *path_;
public:
typedef enum { badEnum,
Open,
High,
Low,
Close,
Volume,
lastEnum } dataKind;
yahooTS()
{
path_ = 0;
};
yahooTS( const char *p ) : path_(p) {}
const double *getTS( const char *fileName,
double *a,
size_t &N,
dataKind kind ) const;
void path( const char *p ) { path_ = p; }
const char *path() { return path_; }
private:
const char *getStr_( char *&line,
char *buf,
size_t bufSize ) const;
void parseVals_( char *line,
double *vals,
const size_t n ) const;
const double getValue_( char *line,
const yahooTS::dataKind kind ) const;
}; // yahooTS
#endif
|
|
|
#include <stdio.h>
#include "yahooTS.h"
main()
{
const size_t VEC_SIZE = 512;
double vec[VEC_SIZE];
yahooTS ts( "equities\\" );
size_t N = VEC_SIZE;
if (ts.getTS( "aa", vec, N, yahooTS::Close )) {
for (size_t i = 0; i < VEC_SIZE; i++) {
printf("%3d %f\n", i, vec[i] );
}
}
else {
fprintf(stderr, "main: getTS failed\n");
}
}
|
|
|
/** \file
The documentation in this file is formatted for doxygen
(see www.doxyeng.org).
<h4>
Copyright and Use
</h4>
<p>
You may use this source code without limitation and without
fee as long as you include:
</p>
<blockquote>
This software was written and is copyrighted by Ian Kaplan, Bear
Products International, www.bearcave.com, 2002.
</blockquote>
<p>
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.
<p>
Please send any bug fixes or suggested source changes to:
<pre>
iank@bearcave.com
</pre>
@author Ian Kaplan
*/
#include <errno.h>
#include <stdio.h>
#include <string.h>
#include "yahooTS.h"
/**
Copy from the input string until either the end of the string
(e.g., the null) is reached or a comma is found.
\param line A reference to a pointer to the input string. This
pointer is incremented until either the end of string or
a comma is encountered. When this function returns
<i>line</i> will either point to the end of the string
or a character following a comma.
\param buf A buffer into which the string will be copied.
\param bufSize The size of <i>buf</i>.
*/
const char *yahooTS::getStr_( char *&line,
char *buf,
unsigned int bufSize ) const
{
const char *rtnPtr = 0;
if (line != 0) {
unsigned int charCnt;
for (charCnt = 0; charCnt < bufSize-1 && *line != '\0'; charCnt++) {
if (*line == ',') {
line++;
break;
}
else {
buf[charCnt] = *line++;
}
}
buf[charCnt] = '\0';
if (charCnt > 0)
{
rtnPtr = buf;
}
}
return rtnPtr;
} // getStr_
/**
Parse a comma separated line of values into a vector of
doubles.
The comma separated values are:
Date,Open,High,Low,Close,Volume
The date value is skipped.
\param line A pointer to a line of Yahoo historical data
\param vals A vector of doubles that the values in the
historical data line will be stored.
\param n The number of elements in <i>vals</i>
*/
void yahooTS::parseVals_( char *line,
double *vals,
const unsigned int n ) const
{
char buf[128];
const char *ptr;
// skip the date
ptr = getStr_( line, buf, sizeof( buf ) );
if (ptr == 0) {
fprintf(stderr, "parseVals: date expected\n" );
return;
}
// get the Open, High, Low, Close and Volume values
unsigned int cnt = 0;
for (dataKind kind = Open;
kind <= Volume && cnt < n;
kind = (dataKind)((unsigned int)kind + 1)) {
ptr = getStr_( line, buf, sizeof( buf ) );
if (ptr == 0) {
fprintf(stderr, "parseVals: value expected\n");
return;
}
double v;
sscanf( buf, "%lf", &v );
vals[cnt] = v;
cnt++;
}
} // parseVals_
/**
A data line from a Yahoo historical data file consists of
a set of comma separated values:
<pre>
date,open,high,low,close,volume
</pre>
This function is passed a Yahoo data line and a <i>kind</i> value
which indicates which value to return. Date is is ignored, so the
value of <i>kind</i> should be one of: Open, High, Low, Close,
Volume.
*/
const double yahooTS::getValue_( char *line,
const yahooTS::dataKind kind ) const
{
double retval = 0;
if (kind > badEnum && kind < lastEnum) {
const unsigned int NUM_VALS = 5;
double vals[ NUM_VALS ];
parseVals_( line, vals, NUM_VALS );
unsigned int ix = (unsigned int)kind - 1;
if (ix < NUM_VALS) {
retval = vals[ix];
}
}
return retval;
} // getValue
/**
Read a Yahoo equity time series from a file.
Yahoo allows historical equity data to be downloaded
in "spread sheet" format. In this format there is a
title line, listing the data columns (e.g., date,
open, high, low, close and volume). Following the
title line are comma separated values. In reading
this Yahoo data file, the first line is skipped.
The Yahoo data values are listed from most recent to
oldest. In the data vector returned, a[0] will be
the oldest and a[N-1] will be the most recent.
\param fileName name of the file containing the time series.
This file will be prefixed by the path in the class
variable path_.
\param a A pointer to a vector of doubles that will
be initialized with values from <i>fileName</i>.
\param N Number of doubles that will fit in <i>a</i>
N is an input/output variable. The value
returned in N will be the actual number of values
read.
\param kind The kind of time series to fetch from
<i>fileName</i> (e.g., open, high, low, close,
volume.
\return If there was no error reading data from <i>fileName</i>
the function returns a pointer to the initialized
array (e.g., the argument <i>a</i>). If the data could
not be read, a null pointer (0) is returned.
*/
const double *yahooTS::getTS( const char *fileName,
double *a,
unsigned int &N,
const yahooTS::dataKind kind ) const
{
const double *rtnPtr = 0;
char fullPath[512];
unsigned int freePath = sizeof( fullPath );
FILE *fptr;
if (path_ != 0) {
strncpy( fullPath, path_, freePath-1 );
freePath = freePath - strlen( fullPath );
}
strncat( fullPath, fileName, freePath-1 );
fptr = fopen( fullPath, "r" );
if (fptr != 0) {
char line[512];
unsigned int lineSize = sizeof( line );
int ix = N-1;
if (fgets( line, lineSize, fptr ) != 0) {
rtnPtr = a;
while (fgets( line, lineSize, fptr ) != 0) {
if (ix >= 0) {
a[ix] = getValue_( line, kind );
ix--;
}
else {
break;
}
} // while
}
else {
fprintf(stderr, "getTS: title line expected\n");
}
ix++;
N = N - ix;
}
else {
const char *error = strerror( errno );
fprintf(stderr, "getTS: Error opening %s: %s\n", fullPath, error );
}
return rtnPtr;
} // getTS
|
|
|
#
# 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
|
|
|
#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).
<h4>
Copyright and Use
</h4>
You may use this source code without limitation and without
fee as long as you include:
<blockquote>
This software was written and is copyrighted by Ian Kaplan, Bear
Products International, www.bearcave.com, 2002.
</blockquote>
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:
<pre>
iank@bearcave.com
</pre>
@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<int> 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 <i>Ripples in
Mathematics</i> 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<int> *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<int> *node) = 0;
public:
/** The default constructor does nothing */
costbase_int() {}
}; // costbase_int
#endif
|
|
|
#ifndef _COSTWIDTH_H_
#define _COSTWIDTH_H_
/** \file
The documentation in this file is formatted for doxygen
(see www.doxygen.org).
<h4>
Copyright and Use
</h4>
You may use this source code without limitation and without
fee as long as you include:
<blockquote>
This software was written and is copyrighted by Ian Kaplan, Bear
Products International, www.bearcave.com, 2002.
</blockquote>
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:
<pre>
iank@bearcave.com
</pre>
@author Ian Kaplan
*/
#include "costbase_int.h"
/**
For a wavelet packet tree node, consisting of <i>n</i> 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<int> *root );
public:
costwidth( packnode<int> *root ) { traverse( root ); }
};
#endif
|
|
|
/**
The delta transform.
A data set of N elements is replaced by the differences. In this
simple algorithm
<pre>
s<sub>i</sub> = s<sub>i</sub> - s<sub>i-1</sub>.
</pre>
For an array indexed from 0, i = 1 to N-1. The element at s<sub>0</sub>
is the reference element, which is unchanged. The next element
at s<sub>1</sub> is replaced by the difference between that element
and s<sub>0</sub>.
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
N<i>log</i><sub>2</sub>N. 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 T>
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, unsigned int len )
{
if (vec != 0 && len > 0) {
// reference value
T next;
T refVal = vec[0];
for (unsigned int i = 1; i < len; i++) {
next = vec[i];
vec[i] = vec[i] - refVal;
refVal = next;
}
}
} // forward
void inverse( T *vec, unsigned int len )
{
if (vec != 0 && len > 0) {
for (unsigned int i = 1; i < len; i++) {
vec[i] = vec[i] + vec[i-1];
}
}
} // inverse
}; // delta_trans
|
|
|
#ifndef _HAAR_INT_H_
#define _HAAR_INT_H_
/** \file
The documentation in this file is formatted for doxygen
(see www.doxygen.org).
<h4>
Copyright and Use
</h4>
<p>
You may use this source code without limitation and without
fee as long as you include:
</p>
<blockquote>
This software was written and is copyrighted by Ian Kaplan, Bear
Products International, www.bearcave.com, 2002.
</blockquote>
<p>
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.
<p>
Please send any bug fixes or suggested source changes to:
<pre>
iank@bearcave.com
</pre>
@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.
<b>References</b>
<ol>
<li>
<p>
<i>Wavelet Transforms that Map Integers to Integers</i> by
A.R. Calderbank, ingrid daubechies, wim weldens and Boon-Lock Yeo,
August 1996
</p>
<p>
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.
</p>
</li>
<li>
<p>
<i>Ripples in Mathematics: the Discrete Wavelet Transform</i>
by Arne Jense and Anders la Cour-Harbo, Springer, 2001
</p>
<p>
This book is a good reference for the Lifting Scheme and the
wavelet transform in general.
</p>
</li>
</ol>
*/
class haar_int : public liftbase<int *, int>
{
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
<a href="http://www.bearcave.com/misl/misl_tech/wavelets/lifting/basiclift.html">Basic Lifting Scheme Wavelets</a>.
*/
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
|
|
|
#ifndef _INVPACKTREE_INT_H_
#define _INVPACKTREE_INt_H_
/** \file
The documentation in this file is formatted for doxygen
(see www.doxygen.org).
<h4>
Copyright and Use
</h4>
<p>
You may use this source code without limitation and without
fee as long as you include:
</p>
<blockquote>
This software was written and is copyrighted by Ian Kaplan, Bear
Products International, www.bearcave.com, 2002.
</blockquote>
<p>
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.
<p>
Please send any bug fixes or suggested source changes to:
<pre>
iank@bearcave.com
</pre>
@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
<i>Ripples in Mathematics</i> 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 <i>Ripples in
Mathematics</i>.
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<packcontainer_int, int> *waveObj;
/** disallow the copy constructor */
invpacktree_int( const invpacktree_int &rhs ) {}
/** inverse wavelet packet transform calculation stack */
LIST<packcontainer_int *> stack;
/** pointer to the inverse transform result */
const int *data;
/** length of data */
unsigned int N;
private:
void new_level( packdata<int> *elem );
void add_elem( packdata<int> *elem );
void reduce();
public:
invpacktree_int( packdata_list<int> &list,
liftbase<packcontainer_int, int> *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
|
|
|
#ifndef _LINE_INT_H_
#define _LINE_INT_H_
/** \file
The documentation in this file is formatted for doxygen
(see www.doxygen.org).
<h4>
Copyright and Use
</h4>
<p>
You may use this source code without limitation and without
fee as long as you include:
</p>
<blockquote>
This software was written and is copyrighted by Ian Kaplan, Bear
Products International, www.bearcave.com, 2002.
</blockquote>
<p>
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.
<p>
Please send any bug fixes or suggested source changes to:
<pre>
iank@bearcave.com
</pre>
@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 T>
class line_int : public liftbase<T, int>
{
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 (b<sub>j</sub>) to the second half
of the array, leaving the even elements (a<sub>i</sub>) in the
first half
<pre>
a<sub>0</sub>, a<sub>1</sub>, a<sub>1</sub>, a<sub>3</sub>, b<sub>0</sub>, b<sub>1</sub>, b<sub>2</sub>, b<sub>2</sub>,
</pre>
The predict step of the line wavelet "predicts" that the
odd element will be on a line between two even elements.
<pre>
b<sub>j+1,i</sub> = b<sub>j,i</sub> - (a<sub>j,i</sub> + a<sub>j,i+1</sub>)/2
</pre>
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 2<sup>n</sup> 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 <tt>new_n_plus1()</tt> 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
(<i>Building Your Own Wavelets at Home</i>) and in <i>Ripples in
Mathematics</i>. A somewhat more complete derivation of this
equation is provided in <i>Ripples in Mathematics</i> 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>i</i>. Note that the predict phase has already
completed, so the odd values belong to iteration <i>i+1</i>.
<pre>
even<sub>i+1,j</sub> = even<sub>i,j</sub> op (odd<sub>i+1,k-1</sub> + odd<sub>i+1,k</sub>)/4
</pre>
This version of the line wavelet code implements an integer
version of linear interpolating wavelet. This versoin comes from
the paper <i>Wavelet Transforms that Map Integers to Integers</i>
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 even<sub>i,j</sub>
(depending on whether the forward or inverse transform is being
calculated) is calculated from odd<sub>i+1,k-1</sub> and
odd<sub>i+1,k</sub> 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
|
|
|
#ifndef _PACKCONTAINER_INT_H_
#define _PACKCONTAINER_INT_H_
#include "packnode.h"
/**
The documentation in this file is formatted for doxygen
(see www.doxygen.org).
<h4>
Copyright and Use
</h4>
<p>
You may use this source code without limitation and without
fee as long as you include:
</p>
<blockquote>
This software was written and is copyrighted by Ian Kaplan, Bear
Products International, www.bearcave.com, 2002.
</blockquote>
<p>
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.
<p>
Please send any bug fixes or suggested source changes to:
<pre>
iank@bearcave.com
</pre>
@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 */
unsigned int 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<int>* node )
{
assert( node != 0 );
N = node->length();
assert( N > 1 );
unsigned int half = N >> 1;
block_pool mem_pool;
unsigned int num_bytes = half * sizeof(int);
lhs = (int *)mem_pool.pool_alloc( num_bytes );
rhs = (int *)mem_pool.pool_alloc( num_bytes );
for (unsigned int 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( unsigned int 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 unsigned int i )
{
assert( i < N );
unsigned int half = N >> 1;
if (i < half)
return lhs[i];
else {
return rhs[i-half];
}
}
/** RHS [] operator */
int operator[]( const unsigned int i ) const
{
assert( i < N );
unsigned int 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.
*/
unsigned int length() { return N; }
}; // packcontainer_int
#endif
|
|
|
#ifndef _PACKTREE_BASE_INT_H_
#define _PACKTREE_BASE_INT_H_
/** \file
The documentation in this file is formatted for doxygen
(see www.doxygen.org).
<h4>
Copyright and Use
</h4>
<p>
You may use this source code without limitation and without
fee as long as you include:
</p>
<blockquote>
This software was written and is copyrighted by Ian Kaplan, Bear
Products International, www.bearcave.com, 2002.
</blockquote>
<p>
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.
<p>
Please send any bug fixes or suggested source changes to:
<pre>
iank@bearcave.com
</pre>
@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<int> *root;
/** wavelet packet transform object */
liftbase<packcontainer_int, int> *waveObj;
typedef enum { BadPrintKind,
printData,
printCost,
printBestBasis } printKind;
void breadthFirstPrint(printKind kind);
void newLevel( packnode<int>* top, bool freqCalc, bool reverse );
public:
void pr();
/** get the root of the wavelet packet tree */
packnode<int> *getRoot() { return root; }
}; // packtree_base_int
#endif
|
|
|
#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).
<h4>
Copyright and Use
</h4>
<p>
You may use this source code without limitation and without
fee as long as you include:
</p>
<blockquote>
This software was written and is copyrighted by Ian Kaplan, Bear
Products International, www.bearcave.com, 2002.
</blockquote>
<p>
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.
<p>
Please send any bug fixes or suggested source changes to:
<pre>
iank@bearcave.com
</pre>
@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
<i>int</i> vaues, the result of the constructor will be a wavelet packet
tree with log<sub>2</sub>(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<int> *root );
void buildBestBasisList( packnode<int> *root,
packdata_list<int> &list );
void checkBestBasis( packnode<int> *root );
void cleanTree(packnode<int> *root, bool removeMark );
public:
packtree_int( const int *vec,
const unsigned int n,
liftbase<packcontainer_int, int> *w );
/** destructor does nothing */
~packtree_int() {}
void prCost();
void prBestBasis();
void bestBasis();
bool bestBasisOK();
packdata_list<int> getBestBasisList();
}; // packtree_int
#endif
|
|
|
#ifndef _SUPPORT_H_
#define _SUPPORT_H_
/** \file
The documentation in this file is formatted for doxygen
(see www.doxygen.org).
<h4>
Copyright and Use
</h4>
<p>
You may use this source code without limitation and without
fee as long as you include:
</p>
<blockquote>
This software was written and is copyrighted by Ian Kaplan, Bear
Products International, www.bearcave.com, 2002.
</blockquote>
<p>
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.
<p>
Please send any bug fixes or suggested source changes to:
<pre>
iank@bearcave.com
</pre>
@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.
<b>Bit width calculation</b>
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.
<b>Double to integer conversion</b>
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 unsigned int nearestPower2Width_( unsigned int 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 unsigned int valWidth( const int val );
unsigned int UnsignedValWidth( const unsigned int val );
static unsigned int vecWidth( const int *vec,
const unsigned int N );
static void roundToInt( int *intVec,
const double *realVec,
const unsigned int len );
static void decimalToInt( int *intVec,
const double *realVec,
const unsigned int len );
}; // support
#endif
|
|
|
#ifndef _TS_TRANS_INT_H_
#define _TS_TRANS_INT_H_
/** \file
The documentation in this file is formatted for doxygen
(see www.doxygen.org).
<h4>
Copyright and Use
</h4>
<p>
You may use this source code without limitation and without
fee as long as you include:
</p>
<blockquote>
This software was written and is copyrighted by Ian Kaplan, Bear
Products International, www.bearcave.com, 2002.
</blockquote>
<p>
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.
<p>
Please send any bug fixes or suggested source changes to:
<pre>
iank@bearcave.com
</pre>
@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.
<pre>
d(1)<sub>1,i</sub> = s<sub>0,2i+1</sub> - s<sub>0,2i</sub>
s<sub>1,i</sub> = s<sub>0,2i</sub> + floor( d(1)<sub>1,i</sub>/2 )
d<sub>1,i</sub> = d(1)<sub>1,i</sub> + floor((s<sub>1,i-1</sub> - s<sub>1,i+1</sub>)/4.0 + 0.5)
</pre>
This notation and the algorithm implemented here is taken directly
from <i>Wavelet Transformst that Map Integers to Integers</i> 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 <tt>predict()</tt> and <tt>update</tt>
functions. An additional <tt>predict2()</tt> function is added
that implements interpolation step.
Since an extra step has been added, the <tt>forwardStep</tt> and
<tt>inverseStep</tt> functions in the base class (<tt>liftbase</tt>)
are over-ridden by ts_trans_int.
<b>A brief note on vanishing moments</b>
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
<pre>
x = {0..2Pi} : y = sin(x)
x = anything else : y = 0
</pre>
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<int *, int>
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 s<sub>1,i+1</sub>.
The low pass half of the array is in the array index range
{0..half-1}. In the equation below, when <tt>i = half-1</tt>
a non-existent element at i+1 is needed. This function
calculates this element from s[half-2] and s[half-1].
<pre>
d<sub>1,i</sub> = d(1)<sub>1,i</sub> + floor((s<sub>1,i-1</sub> - s<sub>1,i+1</sub>)/4.0 + 0.5)
</pre>
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
<pre>
. y<sub>2</sub> - y<sub>1</sub>
(y - y<sub>1</sub>) = -------- (x - x<sub>1</sub>)
. x<sub>2</sub> - x<sub>1</sub>
</pre>
Solving for y
<pre>
. y<sub>2</sub> - y<sub>1</sub>
y = -------- (x - x<sub>1</sub>) + y<sub>1</sub>
. x<sub>2</sub> - x<sub>1</sub>
</pre>
where
<pre>
x<sub>1</sub> = 0
x<sub>2</sub> = 1
y<sub>1</sub> = s[half-2]
y<sub>2</sub> = s[half-1]
x = 2
</pre>
Substituting in these values we get
<pre>
y = 2*y<sub>2</sub> - y<sub>1</sub>
</pre>
*/
int new_n_plus1( int y1, int y2)
{
int y = 2 * y2 - y1;
return y;
}
/**
In the function <tt>new_n_plus1</tt> a point beyond the end
of the low pass filter array is calculated. Here
a point beyond the beginning of the array is calculated.
<pre>
x<sub>1</sub> = 0
x<sub>2</sub> = 1
y<sub>1</sub> = s[0]
y<sub>2</sub> = s[1]
x = -1
</pre>
When these values are plugged into the
two point equation we get
<pre>
y = 2 * y<sub>1</sub> - y<sub>2</sub>
</pre>
*/
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
|
|
|
/** \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).
<h4>
Copyright and Use
</h4>
<p>
You may use this source code without limitation and without
fee as long as you include:
</p>
<blockquote>
This software was written and is copyrighted by Ian Kaplan, Bear
Products International, www.bearcave.com, 2002.
</blockquote>
<p>
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.
<p>
Please send any bug fixes or suggested source changes to:
<pre>
iank@bearcave.com
</pre>
@author Ian Kaplan
*/
#include <stdio.h>
#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 unsigned int N )
{
int *newVec = new int[ N ];
for (unsigned int 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 unsigned int N )
{
bool rslt = true;
for (unsigned int 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.
*/
unsigned int calcPacketWidth( packdata_list<int> &bestBasisList )
{
packdata_list<int>::handle h;
unsigned int totalWidth = 0;
for (h = bestBasisList.first(); h != 0; h = bestBasisList.next( h )) {
packdata<int> *node = bestBasisList.get_item( h );
const unsigned int 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.
*/
unsigned int packet_calc( const int *intVec,
const int *copyVec,
const int N )
{
// The "line" wavelet transform (e.g., line with slope)
line_int<packcontainer_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<int> *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<int> bestBasis = tree.getBestBasisList();
// Sum the cost values (width for each node in the best basis
// list
unsigned int 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 s<sub>i-1</sub> and s<sub>i</sub>.
*/
unsigned int delta_calc( int *intVec,
const int *copyVec,
const int N )
{
delta_trans<int> delta;
delta.forward( intVec, N );
const unsigned int 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
*/
unsigned int wave_calc( int *intVec,
const int *copyVec,
const int N,
liftbase<int *, int> *w )
{
w->forwardTrans( intVec, N );
const unsigned int 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 unsigned int N = 512;
double realVec[ N ];
int intVec[ N ];
// an instance of yahooTS with the path to the data directory
const char *dataDirPath = "..\\haardata\\equities\\";
yahooTS ts( dataDirPath );
printf("Equity Uncompressed delta Haar line TS wavelet packet (line)\n");
for (unsigned int i = 0; files[i] != 0; i++) {
unsigned int 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 unsigned int beforeWidth = support::vecWidth( intVec, N );
const unsigned int deltaWidth = delta_calc( intVec, copyVec, N );
haar_int haar;
const unsigned int haarWidth = wave_calc( intVec, copyVec, N, &haar );
line_int<int *> line;
const unsigned int lineWidth = wave_calc( intVec, copyVec, N, &line );
ts_trans_int ts_trans;
const unsigned int tsTransWidth = wave_calc( intVec,
copyVec,
N,
&ts_trans);
const unsigned int 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;
}
|
|
|
/** \file
The documentation in this file is formatted for doxygen
(see www.doxygen.org).
<h4>
Copyright and Use
</h4>
You may use this source code without limitation and without
fee as long as you include:
<blockquote>
This software was written and is copyrighted by Ian Kaplan, Bear
Products International, www.bearcave.com, 2002.
</blockquote>
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:
<pre>
iank@bearcave.com
</pre>
@author Ian Kaplan
*/
#include "costwidth.h"
#include "support.h"
int
costwidth::costCalc( packnode<int> *root )
{
assert( root != 0 );
unsigned int N = root->length();
const int *a = root->getData();
int width = support::vecWidth( a, N );
return width;
} // costCalc
|
|
|
#include <assert.h>
#include <stdio.h>
#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).
<h4>
Copyright and Use
</h4>
<p>
You may use this source code without limitation and without
fee as long as you include:
</p>
<blockquote>
This software was written and is copyrighted by Ian Kaplan, Bear
Products International, www.bearcave.com, 2002.
</blockquote>
<p>
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.
<p>
Please send any bug fixes or suggested source changes to:
<pre>
iank@bearcave.com
</pre>
@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<int> *elem )
{
unsigned int half = elem->length();
unsigned int 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<packcontainer_int *>::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();
unsigned int 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<int> *elem )
{
assert( elem != 0 );
if (stack.first() == 0) {
new_level( elem );
}
else {
unsigned int half = elem->length();
unsigned int n = half * 2;
LIST<packcontainer_int *>::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<int> &list,
liftbase<packcontainer_int, int> *w )
{
data = 0;
N = 0;
waveObj = w;
// Traverse the "best basis" list and calculate the inverse
// wavelet packet transform.
packdata_list<int>::handle h;
for (h = list.first(); h != 0; h = list.next( h )) {
packdata<int> *elem = list.get_item( h );
add_elem( elem );
} // for
LIST<packcontainer_int *>::handle tosHandle;
tosHandle = stack.first();
packcontainer_int *tos = stack.get_item( tosHandle );
if (tos != 0) {
unsigned int 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
|
|
|
#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 <i>top</i> packnode object contains data from the previous
level. If this is the first level in the tree, <i>top</i> will
contain the input data set.
A packcontainer_int object is created with the <i>top</i> 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 <i>top</i>. 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<int>* top,
bool freqCalc,
bool reverse )
{
if (top != 0) {
const unsigned int 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<int> *lhs = new packnode<int>(container.lhsData(),
len/2,
packnode<int>::LowPass );
packnode<int> *rhs = new packnode<int>(container.rhsData(),
len/2,
packnode<int>::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 <i>Fundamentals of Data Structures in C</i> by Horowitz, Sahni
and Anderson-Freed, 1993.
*/
void packtree_base_int::breadthFirstPrint(const printKind kind)
{
queue<int> Q;
Q.addQueue( root, 0 );
while (! Q.queueEmpty() ) {
packnode<int> *node = Q.queueStart()->node;
unsigned int 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<int> *lhs = node->lhsChild();
packnode<int> *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
|
|
|
/** \file
The documentation in this file is formatted for doxygen
(see www.doxyeng.org).
<h4>
Copyright and Use
</h4>
<p>
You may use this source code without limitation and without
fee as long as you include:
</p>
<blockquote>
This software was written and is copyrighted by Ian Kaplan, Bear
Products International, www.bearcave.com, 2002.
</blockquote>
<p>
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.
<p>
Please send any bug fixes or suggested source changes to:
<pre>
iank@bearcave.com
</pre>
@author Ian Kaplan
*/
#include <assert.h>
#include <stdio.h>
#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 <i>w</i>. 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 unsigned int N,
liftbase<packcontainer_int, int> *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<int>( vecCopy, N, packnode<int>::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 <i>Ripples in Mathematics</i>
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<int> *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 <i>Ripples in Mathematics</i>. 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<int> *top )
{
int cost = 0;
if (top != 0) {
packnode<int> *lhs = top->lhsChild();
packnode<int> *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<int> *top )
{
if (top != 0) {
if (top->mark()) {
foundBestBasisVal = true;
if (top->getKind() == packdata<int>::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<int> *top,
packdata_list<int> &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<int> packtree_int::getBestBasisList()
{
packdata_list<int> list;
buildBestBasisList( root, list );
return list;
} // getBestBasisList
|
|
|
#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:
<ul>
<li>
If the digit in the fourth
place is greater than 5, the digit in the third place is incremented
by one.
</li>
<li>
If the digit in the fourth place is less than 5, the third digit
is left unchanged.
</li>
<li>
If the digit in the fourth place is 5 and the digit in the third
place is odd, the third digit is incremented.
</li>
</ul>
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:
<pre>
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.
</pre>
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<int>(val);
double fracPart = val - intPart;
int threeDigits = static_cast<int>(fracPart * 1000.0);
int fourDigits = static_cast<int>(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 static_cast<int>(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 unsigned int len )
{
if (intVec != 0 && realVec != 0) {
for (unsigned int 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 <tt>roundToInt</tt> function can be used for rounding.
*/
void support::decimalToInt( int *intVec,
const double *realVec,
const unsigned int len )
{
if (intVec != 0 && realVec != 0) {
for (unsigned int i = 0; i < len; i++) {
intVec[i] = (int)(realVec[i] * 100.0);
}
}
} // roundToInt
unsigned int support::nearestPower2Width_( unsigned int val )
{
unsigned int width = 0;
if (val > 0) {
width = 1;
unsigned int 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.
*/
unsigned int support::valWidth( const int val )
{
unsigned int wholeNum = (val < 0) ? -val : val;
unsigned int width = 1 + nearestPower2Width_( wholeNum );
return width;
} // valWidth
/**
Calculate the number of bits needed to represent an unsigned
value.
*/
unsigned int support::UnsignedValWidth( const unsigned int val )
{
unsigned int width = nearestPower2Width_( val );
return width;
} // valWidth
/**
Calculate the minimum number of bits needed to represent the
values in an integer vector.
*/
unsigned int support::vecWidth( const int *vec, const unsigned int N )
{
unsigned int totalWidth = 0;
if (vec != 0) {
for (unsigned int i = 0; i < N; i++) {
totalWidth += valWidth( vec[i] );
}
}
return totalWidth;
} // vecWidth( int *)
|