The Daubechies D4 Wavelet Transform in C++ and Java

http://www.bearcave.com/software/java/wavelets/daubechies/index.html

I recommend using the "save as" feature of your browser to save the C++ and Java source files (I'm not sure how to reliably suppress viewing and force download with all browsers).

• The Daubechies D4 algorithm as a C++ class

• Doxygen generated documentation for the C++ version of the Daubechies D4 algorithm. This documentation was generated with this Doxygen configuration file.

• The Daubechies D4 algorithm as a Java class

Daubechies D4 wavelet transform (D4 denotes four coefficients)

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

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

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

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

depending on the version of the Haar algorithm used.

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

This comment is formatted for the doxygen documentation generator

Reading and writing wav files

This example demonstrates how to use Dr. Fred DePiero's classes to read and write to/from files in the .wav format. All of Dr. DePiero's classes needed to read/write .wav files are included without comment. The program also demonstrates the use of command line arguments.

download: http://yongzhen.zhuang.googlepages.com/readwritewavfiles.rar

Pseudo random number generators - uniform and non-uniform distributions

perfect one!!!! contain many generators!!!

http://www.agner.org/random/

This page contains software libraries for some very good random number generators.

The basic random number generators make floating point or integer random numbers with uniform distributions. This code is available in C++ and assembly language.

The non-uniform random number generators make random variates with the following distributions:
normal, bernoulli, poisson, binomial, hypergeometric, Wallenius' and Fisher's noncentral hypergeometric, multinomial, multivariate hypergeometric, Wallenius' and Fisher's multivariate noncentral hypergeometric, and shuffling. This code is available in C++ language.

Code examples are included, showing how to use these software libraries.

The uniform random number generators are also available as ready-to-use library files which can be linked into projects in many different programming languages under Windows, Linux, BSD and other operating systems on the PC platform. These libraries are coded in assembly language for improved speed.

These generators are intended for Monte Carlo applications, not for cryptographic applications.

Download packages:

Uniform random number generators in C++

Description: C++ class library containing the following random number generators: Mersenne twister and Mother-of-all. Can generate floating point or integer random numbers with uniform distribution, and random bits. Very good randomness, high resolution, extremely long cycle lengths, and high speed. Example included.
Also includes assembly language implementations for x86-based systems for improved speed.
System requirements: Any C++ compiler, any operating system.
Further description and instructions
Description of Mother-of-all generator
File name: randomc.zip, size: 102418, last modified: 2007-Sep-23.
Download C++ random number generators.

Non-uniform random number generators in C++

Description: C++ class library generating random numbers with the following distributions: normal, bernoulli, poisson, binomial, hypergeometric, Wallenius' and Fisher's noncentral hypergeometric, multinomial, multivariate hypergeometric, and multivariate Fisher's and Wallenius' noncentral hypergeometric distributions. A function for shuffling numbers is also included, as well as C++ examples showing how to use these functions for simulating evolution and for other purposes.
Most of the functions are fast and accurate, even for extreme values of the parameters.
You have the choice of using any of the uniform random number generators listed above (C++ or assembly) as base for these non-uniform random number generators.
Further description and instructions.

Definition of distributions pdf format.
Wallenius' noncentral hypergeometric distribution theory.
Theoretical description of sampling methods used pdf format.
File name: stocc.zip, size: 262823, last modified: 2007-Sep-23.
Download non-uniform random number generators.

List of random numbers

Description: A list of 10000 random numbers generated with a combined generator.
File name: 10000ran.zip, size: 49098, last modified: 2005-May-24.
Download 10000 random numbers.

R package for noncentral hypergeometric distributions

Description: Package for the R language (www.r-project.org) for calculating the various noncentral hypergeometric distributions. Useful for biased urn models, models of biased sampling and evolution by natural selection.
Package name: BiasedUrn, last modified: 2007-Jun-16.
BiasedUrn.

Comments to the theory of these random number generators can be posted to my discussion board.

Follow my research on the noncentral hypergeometric distributions.

Please don't mail me with your programming problems. Your mail will not be answered.

Links to related sites

• PLAB: A Server on the Theory and Practice of Random Number Generation
  
• The R Project for Statistical Computing. Free software system for mathematical and statistical computing.
  
• Ernst Stadlober's page on non-uniform random number generation, including the Win-rand package which contains many other distributions than the ones supplied here, and a graphic interface to show the distribution curves.
  
• Ranlib: library of functions for random number generators and non-uniform variates ftp://odin.mda.uth.tmc.edu/pub/source and ftp://odin.mda.uth.tmc.edu/pub/msdos
  
• RngPack: Random number generators for Java: http://www.honeylocust.com/RngPack/
  
• Matsumoto's Mersenne Twister RNG: http://www.math.sci.hiroshima-u.ac.jp/~m-mat/eindex.html
  
• Regress: Statistics package for Macintosh: http://www.causaScientia.org/software/Regress_plus.html
  

use 'printf' in matlab

http://yongzhen.zhuang.googlepages.com/printfinmatlab.rar

random number generator

1. sources

• C++ RNG
• Matlab RNG
• http://www-personal.umich.edu/~wagnerr/MersenneTwister.html

2. tutorial

Intro

This tutorial provides a brief introduction to the random number functions that come as part of the C++ standard library, namely rand() and srand().

rand() and RAND_MAX

The C++ standard library includes a pseudo random number generator for generating random numbers. In order to use it we need to include the header. To generate a random number we use the rand() function. This will produce a result in the range 0 to RAND_MAX, where RAND_MAX is a constant defined by the implementation.

Here's a piece of code that will generate a single random number:

#include 
#include 

using namespace std;

int main()
{
 int random_integer = rand();
 cout <<>

The value of RAND_MAX varies between compilers and can be as low as 32767, which would give a range from 0 to 32767 for rand(). To find out the value of RAND_MAX for your compiler run the following small piece of code:

#include 
#include 

using namespace std;

int main()
{
 cout << "The value of RAND_MAX is " <<>

srand()

The pseudo random number generator produces a sequence of numbers that gives the appearance of being random, when in fact the sequence will eventually repeat and is predictable.

We can seed the generator with the srand() function. This will start the generator from a point in the sequence that is dependent on the value we pass as an argument. If we seed the generator once with a variable value, for instance the system time, before our first call of rand() we can generate numbers that are random enough for simple use (though not for serious statistical purposes).

In our earlier example the program would have generated the same number each time we ran it because the generator would have been seeded with the same default value each time. The following code will seed the generator with the system time then output a single random number, which should be different each time we run the program.

#include 
#include 
#include 

using namespace std;

int main()
{
 srand((unsigned)time(0)); 
 int random_integer = rand();
 cout <<>

Don't make the mistake of calling srand() every time you generate a random number; we only usually need to call srand() once, prior to the first call to rand().

Generating a number in a specific range

If we want to produce numbers in a specific range, rather than between 0 and RAND_MAX, we can use the modulo operator. It's not the best way to generate a range but it's the simplest. If we use rand()%n we generate a number from 0 to n-1. By adding an offset to the result we can produce a range that is not zero based. The following code will produce 20 random numbers from 1 to 10:

#include 
#include 
#include 

using namespace std;

int main()
{
 srand((unsigned)time(0));
 int random_integer;
 for(int index=0; index<20; random_integer =" (rand()%10)+1;">

A better method, though slightly more complicated, is given below. This overcomes problems that are sometimes experienced with some types of pseudo random number generator that might be supplied with your compiler. As before, this will output 20 random numbers from 1 to 10.

#include 
#include 
#include 

using namespace std;

int main()
{
 srand((unsigned)time(0));
 int random_integer;
 int lowest=1, highest=10;
 int range=(highest-lowest)+1;
 for(int index=0; index<20; random_integer =" lowest+int(range*rand()/(RAND_MAX">

Conclusion

If you need to use a pseudo random number generator for anything even remotely serious you should avoid the simple generator that comes with your compiler and use something more sophisticated instead. That said, rand() still has its place and you may find it useful.

kdtree by Steven Michael

Name: kdtree
Author: Steven Michael (smichael@ll.mit.edu)
Date: 3/1/2005

############################################################

The following code implements a KDTree search algorithm
in MATLAB

There are 4 main functions:

1. kdtree -- tree class creation
2. kdtree_range -- return all points within a range
3. kdtree_closestpoint -- return array of closest points to a
corresponding array of input points


A single reference was used in writing the code:

M. deBerg, M. vanKreveld, M. Overmars, and O. Schwarzkopf.
"Computational Geometry: Algorithms and Applications"
Springer, 2000.

download

Anyoption - C/C++ Command line and resource file option parsing

http://www.hackorama.com/anyoption/

AnyOption is a C++ class for easy parsing of complex commandline options. It also parses options from a rsourcefile in option value pair format.

AnyOption implements the traditional POSIX style character options ( -n ) as well as the newer GNU style long options ( --name ). Or you can use a simpler long option version ( -name ) by asking to ignore the POSIX style options.

AnyOption supports the traditional UNIX resourcefile syntax of, any line starting with "#" is a comment and the value pairs use ":" as a delimiter.

An option which expects a value is considered as an option value pair, while options without a value are considered flags.

Please read the header file for the documented public interface, and demo.cpp for an example of how easy it is to use AnyOption.

August 2004, added bug-fixes, and updates send by Michael Peters of Sandia Lab.
September 2006, fix from Boyan Asenov for a bug in mixing up option type indexes.

View Code:
anyoption.cpp
anyoption.h
demo.cpp

Get Code:
anyoption.cpp
anyoption.h
demo.cpp