lifttest.java 100644 1040 1001 4075 7427050360 11763 0 ustar everyone
import lift.*;
import util.*;
class lifttest {
public static void main( String[] args ) {
double vals[] = { 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 vals[] = {
77.6875, 78.1875, 82.0625, 85.5625, 86.7500, 82.4375,
82.2500, 82.7500, 81.2500, 79.5625, 80.2813, 79.8750,
77.7500, 74.7500, 78.5000, 79.1875, 78.8125, 80.3125,
80.1250, 79.3125, 83.7500, 89.8125, 87.7500, 91.1250,
94.4375, 92.7500, 98.0000, 97.1875, 99.4375, 101.7500,
108.5000, 109.0000, 105.2500, 105.5000, 110.0000, 107.0000,
107.2500, 103.3125, 102.8750, 102.4375, 102.0000, 101.3125,
97.4375, 100.5000, 107.7500, 110.2500, 114.3125, 111.2500,
114.8125, 112.6875, 109.4375, 108.0625, 104.5625, 103.2500,
110.5625, 110.7500, 116.3125, 123.6250, 120.9375, 121.6250,
127.6875, 126.0625, 126.3750, 124.3750
};
*/
print pr = new print();
haar hr = new haar();
line ln = new line();
haarpoly hrpy = new haarpoly();
poly py = new poly();
System.out.println("Data:");
pr.pr_vec( vals );
System.out.println();
System.out.println("Haar:");
hr.forwardTrans( vals );
pr.pr_ordered( vals );
System.out.println();
hr.inverseTrans( vals );
pr.pr_vec( vals );
System.out.println();
System.out.println("Line:");
ln.forwardTrans( vals );
pr.pr_ordered( vals );
System.out.println();
ln.inverseTrans( vals );
pr.pr_vec( vals );
System.out.println();
System.out.println("Haar, extended with polynomial interpolation:");
hrpy.forwardTrans( vals );
pr.pr_ordered( vals );
System.out.println();
hrpy.inverseTrans( vals );
pr.pr_vec( vals );
System.out.println();
System.out.println("Poly:");
py.forwardTrans( vals );
pr.pr_ordered( vals );
System.out.println();
py.inverseTrans( vals );
pr.pr_vec( vals );
System.out.println();
}
}
lifttest.cpp 100644 1040 1001 1226 7426143205 11617 0 ustar everyone
#include
The polynomial interpolation algorithm assumes that the known
points are located at x-coordinates 0, 1,.. N-1. An interpolated
point is calculated at x, using N coefficients. The
polynomial coefficients for the point x can be
calculated staticly, using the Lagrange method.
For a given N-point polynomial interpolation, fill the coefficient
table, for points 0.5 ... (N-0.5).
For the polynomial interpolation point x-coordinate
x, return the associated polynomial
interpolation coefficients.
Given four points at the x,y coordinates {0,d0},
{1,d1}, {2,d2}, {3,d3}
return the y-coordinate value for the polynomial interpolated
point at x.
Copy four points or N (which ever is less) data points from
vec into d These points are the "known" points used
in the polynomial interpolation.
Predict an odd point from the even points, using 4-point
polynomial interpolation.
The four points used in the polynomial interpolation are
the even points. We pretend that these four points
are located at the x-coordinates 0,1,2,3. The first
odd point interpolated will be located between the first
and second even point, at 0.5. The next N-3 points are
located at 1.5 (in the middle of the four points).
The last two points are located at 2.5 and 3.5.
For complete documentation see
The difference between the predicted (interpolated) value
and the actual odd value replaces the odd value in the
forward transform.
As the recursive steps proceed, N will eventually be 4 and
then 2. When N = 4, linear interpolation is used.
When N = 2, Haar interpolation is used (the prediction for
the odd value is that it is equal to the even value).
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.
Haar (flat line) wavelet.
As with all Lifting scheme wavelet transform functions, the
first stage of a transform step is the split stage. The
split step moves the even element to the first half of an
N element region and the odd elements to the second half of the N
element region.
The Lifting Scheme version of the Haar transform uses a wavelet
function (predict stage) that "predicts" that an odd element will
have the same value as it preceeding even element. Stated another
way, the odd element is "predicted" to be on a flat (zero slope
line) shared with the even point. The difference between this
"prediction" and the actual odd value replaces the odd element.
The wavelet scaling function (a.k.a. smoothing function) used
in the update stage calculates the average between an even and
an odd element.
The merge stage at the end of the inverse transform interleaves
odd and even elements from the two halves of the array
(e.g., ordering them even0, odd0,
even1, odd1, ...)
You may use this source code without limitation and without
fee as long as you include:
This software is provided "as is", without any warrenty 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:
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
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
Haar transform extended with a polynomial interpolation step
This wavelet transform extends the Haar transform with a
polynomial wavelet function.
The polynomial wavelet uses 4-point polynomial interpolation to
"predict" an odd point from four even point values.
This class extends the Haar transform with an interpolation
stage which follows the predict and update stages of the
Haar transform. The predict value is calculated from the
even points, which in this case are the smoothed values
calculated by the scaling function (e.g., the averages
of the even and odd values).
The predict value is subtracted from the current odd value, which
is the result of the Haar wavelet function (e.g., the difference
between the odd value and the even value). This tends to
result in large odd values after the interpolation stage, which
is a weakness in this algorithm.
This algorithm was suggested by Wim Sweldens' tutorial
Building Your Own Wavelets at Home.
f(x) = c0y0 + c1y1 + c2y2 + c2y2;
Where f(x) is the interpolated point and
{y0, y1, y2, y3} are
four points from a time series.
The reference I used for this is Chapter 3.1 of Numerical
Recipes in C. In Numerical Recipies they criticize the Lagrange
polynomial method because it does not calculate an error value
for the calculated point. For some sets of four points the
value calculated will have a high error and is unreliable. But
for the purposes of wavelets the error is not really useful.
Calculating the coefficients in advance for a four point
Lagrange polynomial at "x" is faster than using Neville's
algorithm every time.
*/
void lagrange( double coef[], double x )
{
double x1 = 0;
double x2 = 1;
double x3 = 2;
double x4 = 3;
coef[0] = ((x - x2)* (x - x3) * (x - x4)) / ((x1 - x2) * (x1 - x3) * (x1 - x4));
coef[1] = ((x - x1)* (x - x3) * (x - x4)) / ((x2 - x1) * (x2 - x3) * (x2 - x4));
coef[2] = ((x - x1)* (x - x2) * (x - x4)) / ((x3 - x1) * (x3 - x2) * (x3 - x4));
coef[3] = ((x - x1)* (x - x2) * (x - x3)) / ((x4 - x1) * (x4 - x2) * (x4 - x3));
} // lagrange
/**
Calculate the polynomial interpolation for a point x.
The function is passed four data points in vals and
a four value coefficient array, containing the coefficients for
the polymomial at point x (for example 1.5). This function
calculates
f(x) = c0y0 + c1y1 + c2y2 + c2y2;
*/
double lagrangeCalc( double vals[], double coef[] )
{
double rslt = (vals[0] * coef[0]) +
(vals[1] * coef[1]) +
(vals[2] * coef[2]) +
(vals[3] * coef[3]);
return rslt;
} // lagrangeCalc
/**
Print the polynomial coefficients. This is used for debugging.
*/
void pr_coef( double coef[] )
{
printf("{%7.4f, %7.4f, %7.4f, %7.4f}\n",
coef[0], coef[1], coef[2], coef[3] );
}
public:
/**
The interpolation class constructor calculates the coefficients for
a four point polynomial interpolated at -0.5, 0.5, 1.5, 2.5 and 3.5.
*/
interpolation()
{
lagrange( coef_minus2, -2.0 );
lagrange( coef_minus1, -1.0 );
lagrange( coef_minus05, -0.5 );
lagrange( coef_05, 0.5 );
lagrange( coef_15, 1.5 );
lagrange( coef_25, 2.5 );
lagrange( coef_35, 3.5 );
lagrange( coef_40, 4.0 );
lagrange( coef_50, 5.0 );
} // interpolation class constructor
/**
For a point x calculate the polynomial interpolation
from the set of points in vals. The argument n
is the size of the time series (averages). Note that as
the wavelet transform calculation recursively proceeds, n
decreases by a factor of two each time (e.g., 32, 16, 8...).
*/
double pointCalc( double vals[], const double x, const int n )
{
double rslt = 0.0;
if (n == 2) {
if (x == -0.5) {
rslt = (vals[0] * (-0.5)) + (vals[1] * 1.5);
}
else if (x == 0.5) {
// interval betwee 0 .. 1
rslt = (vals[0] * 0.5) + (vals[1] * 0.5);
}
else if (x == 1.5) {
// interval between 1 .. 2
rslt = (vals[0] * 1.5) + (vals[1] * (-0.5));
}
else
printf("pointCalc, n = 2: can't calculate point at that x = %7.4f\n", x);
}
else { // n >= 4
int ix = (int)(x * 10);
switch (ix) {
case -20: // -2.0
{
rslt = lagrangeCalc( vals, coef_minus2 );
}
break;
case -10: // -1.0
{
rslt = lagrangeCalc( vals, coef_minus1 );
}
break;
case -5: // -0.5
{
rslt = lagrangeCalc( vals, coef_minus05 );
}
break;
case 5: // 0.5
{
rslt = lagrangeCalc( vals, coef_05 );
}
break;
case 15: // 1.5
{
rslt = lagrangeCalc( vals, coef_15 );
}
break;
case 25: // 2.5
{
rslt = lagrangeCalc( vals, coef_25 );
}
break;
case 35: // 3.5
{
rslt = lagrangeCalc( vals, coef_35 );
}
break;
case 40: // 4.0
{
rslt = lagrangeCalc( vals, coef_40 );
}
break;
case 50: // 5.0
{
rslt = lagrangeCalc( vals, coef_50 );
}
break;
default:
{
printf("pointCalc, n >= 4: can't calculate point at that x = %7.4f\n", x);
}
break;
} // switch
}
return rslt;
} // pointCalc
}; // class interpolation
C++/leastsq.h 100644 1040 1001 3055 7363726044 11424 0 ustar everyone
/**
Least squares fit for four evenly spaced points on the y-axis.
The x-axis values are, respectively, {0, 1, 2, 3}. As with
polynomial interpolation for four points at known x-axis
locations degenerates into a set of linear coefficients.
x0..x3 = {0, 1, 2, 3}. So the mean of the x-axis points is
Xm = 6/4 = 1.5.
The sum of the squares around the mean for X is similarly
constant since we know that Xi = {0, 1, 2, 3} for the
four points.
SSxx = sq(0 - 1.5) + sq(1 - 1.5) + sq(2 - 1.5) + sq(3 - 1.5)
SSxx = sq(-1.5) + sq(-0.5) + sq(0.5) + sq(1.5)
SSxx = 2.25 + 0.25 + 0.25 + 2.25
SSxx = 5.0
For the linear equation
y = a + bx
b = SSxy / SSxx
Ym = mean of Yi
SSxy = (-1.5)(y0 - Ym) + (-0.5)(y1 - Ym) + (0.5)(y2 - Ym) + (1.5)(y3 - Ym)
SSxy = -3(y0 - Ym)/2 + -(y1 - Ym)/2 + (y2 - Ym)/2 + 3(y3 - Ym)/2
SSxy = (-3y0 + 3Ym)/2 + (-y1 + Ym)/2 + (y2 - Ym)/2 + (3y3 - 3Ym)/2
SSxy = (-3y0 + 3Ym + (-y1) + Ym + y2 - Ym + 3y3 - 3Ym)/2
SSxy = (-3y0 - y1 + y2 + 3y3 )/2
b = SSxy / SSxx = ((-3y0 - y1 + y2 + 3y3 )/2)/5
b = SSxy / SSxx = (-3y0 - y1 + y2 + 3y3 )/10
*/
class leastsq {
private:
double a, b;
public:
leastsq( double y0, double y1, double y2, double y3 )
{
const double Xm = 1.5;
double Ym = (y0 + y1 + y2 + y3)/4.0;
b = ((-3 * y0) - y1 + y2 + (3 * y3))/10.0;
a = Ym - (b * Xm);
}
double yPoint( double x )
{
double y = a + b * x;
return y;
}
};
C++/lifting.h 100644 1040 1001 43246 7426422112 11417 0 ustar everyone
#include
http://www.bearcave.com/misl/misl_tech/wavelets/lifting/index.html
a0 b0 a1 b1 a2 b2 a3 b3
The split function will move the odd elements into the second half
of the array, resulting in
a0 a2 a3 a4 b0 b2 b3 b4
*/
void split( DoubleVec& ts, const int N )
{
int start = 1;
int end = N-1;
while (start < end) {
int i;
double tmp;
for (i = start; i < end; i = i + 2) {
tmp = ts[i];
ts[i] = ts[i+1];
ts[i+1] = tmp;
}
start = start + 1;
end = end - 1;
}
} // split
/**
Merge the odd and even values. The is the inverse of the
split.
If the even values are labeled with a and the odd values
are labeled with b, the input vector will contain
a0 a2 a3 a4 b0 b2 b3 b4
The merge function will re-order the odd and even element in the
vector:
a0 b0 a1 b1 a2 b2 a3 b3
*/
void merge( DoubleVec& ts, int N)
{
int half = N >> 1;
int start = half-1;
int end = half;
int i;
double tmp;
while (start > 0) {
for (i = start; i < end; i = i + 2) {
tmp = ts[i];
ts[i] = ts[i+1];
ts[i+1] = tmp;
}
start = start - 1;
end = end + 1;
}
} // merge
/**
Predict step of the wavelet Lifting Scheme.
The term "predict" comes from Building Your Own Wavelets
at Home.
transform
In the wavelet transform, calculate the difference between even
and odd element pairs. This is a slightly modified version of
the classic haar difference. If the even elements are labeled
as a and the odd elements are labeled as b (where
we start counting at zero) the difference is simply
di = bi - ai
Since an "in-place" algorithm is used, where we update the
odd elements with the difference, we get
bi = bi - ai
inverse transform
Reverse the transform predict step by adding the average
(even element) to the difference (odd element).
*/
void predict( DoubleVec& ts, const int N, const direction dir)
{
int i;
int half = N >> 1;
int j = 0;
for (i = half; i < N; i++) {
if (dir == forward) { // forward transform stage
ts[i] = ts[i] - ts[j];
}
else if (dir == inverse ) { // inverse transform stage
ts[i] = ts[i] + ts[j];
}
else {
printf("predict: bad direction\n");
break;
}
j++;
}
} // predict
/**
The update step of the wavelet Lifting Scheme
transform
In the Lifting Scheme transform the update step follows
the predict step. After the predict step, the differences
have been calculated in-place writing over the even (b)
values. The update step calculates the Haar average using
an even/odd pair. The average is placed in the even member
of the pair.
*/
void update( DoubleVec& ts, int N, const direction dir)
{
int i;
int half = N >> 1;
int j = half;
for (i = 0; i < half; i++) {
if (dir == forward) { // forward transform stage
ts[i] = ts[i] + (ts[j]/2.0);
}
else if (dir == inverse) { // inverse transform step
ts[i] = ts[i] - (ts[j]/2.0);
}
else {
printf("update: bad direction value\n");
break;
}
j++;
} // for
} // update
public:
/**
Wavelet lifting scheme transform.
The wavelet lifting scheme is an in-place wavelet algorithm.
A time series of N elements is passed to the transform
function. The result of the wavelet lifting scheme is calculated
in place (without any array temporaries) in the ts
DoubleVec.
*/
void transform( DoubleVec& ts )
{
const int N = ts.size();
int n;
for (n = N; n > 1; n = n >> 1) {
split( ts, n );
predict( ts, n, forward );
update( ts, n, forward );
interp( ts, n, forward );
} // for
} // transform
/**
Wavelet lifting Scheme inverse transform.
Like the Lifting Scheme transform, the inverse transform is an
in-place algorithm. The result of the transform is passed to
the inverse transform, which calculates the time series from
the time series average and the coefficients.
*/
void invTransform( DoubleVec& ts )
{
const int N = ts.size();
int n;
for (n = 2; n <= N; n = n << 1) {
interp( ts, n, inverse );
update( ts, n, inverse );
predict( ts, n, inverse );
merge( ts, n );
}
} // invTransform
}; // class liftingScheme
lift/ 40755 1040 1001 0 7427050163 10116 5 ustar everyone lift/daubechies.java 100644 1040 1001 6761 7424663416 13174 0 ustar everyone
package lift;
import java.lang.Math.*;
public class daubechies extends liftbase {
final static double sqrt3 = Math.sqrt( 3 );
final static double sqrt2 = Math.sqrt( 2 );
protected void normalize( double[] S, int N, int direction )
{
int half = N >> 1;
for (int n = 0; n < half; n++) {
if (direction == forward) {
S[n] = ((sqrt3-1.0)/sqrt2) * S[n];
S[n+half] = ((sqrt3+1.0)/sqrt2) * S[n+half];
}
else if (direction == inverse) {
S[n] = ((sqrt3+1.0)/sqrt2) * S[n];
S[n+half] = ((sqrt3-1.0)/sqrt2) * S[n+half];
}
else {
System.out.println("daubechies::normalize: bad direction value");
break;
}
}
} // normalize
protected void predict( double[] S, int N, int direction )
{
int half = N >> 1;
if (direction == forward) {
S[half] = S[half] - (sqrt3/4.0)*S[0] - (((sqrt3-2)/4.0)*S[half-1]);
}
else if (direction == inverse) {
S[half] = S[half] + (sqrt3/4.0)*S[0] + (((sqrt3-2)/4.0)*S[half-1]);
}
else {
System.out.println("daubechies::predict: bad direction value");
}
// predict, forward
for (int n = 1; n < half; n++) {
if (direction == forward) {
S[half+n] = S[half+n] - (sqrt3/4.0)*S[n] - (((sqrt3-2)/4.0)*S[n-1]);
}
else if (direction == inverse) {
S[half+n] = S[half+n] + (sqrt3/4.0)*S[n] + (((sqrt3-2)/4.0)*S[n-1]);
}
else {
break;
}
}
} // predict
protected void updateOne( double[] S, int N, int direction )
{
int half = N >> 1;
for (int n = 0; n < half; n++) {
double updateVal = sqrt3 * S[half+n];
if (direction == forward) {
S[n] = S[n] + updateVal;
}
else if (direction == inverse) {
S[n] = S[n] - updateVal;
}
else {
System.out.println("daubechies::updateOne: bad direction value");
break;
}
}
} // updateOne
protected void update( double[] S, int N, int direction )
{
int half = N >> 1;
for (int n = 0; n < half-1; n++) {
if (direction == forward) {
S[n] = S[n] - S[half+n+1];
}
else if (direction == inverse) {
S[n] = S[n] + S[half+n+1];
}
else {
System.out.println("daubechies::update: bad direction value");
break;
}
}
if (direction == forward) {
S[half-1] = S[half-1] - S[half];
}
else if (direction == inverse) {
S[half-1] = S[half-1] + S[half];
}
} // update
public void forwardTrans( double[] vec )
{
final int N = vec.length;
for (int n = N; n > 1; n = n >> 1) {
split( vec, n );
updateOne( vec, n, forward ); // update 1
predict( vec, n, forward );
update( vec, n, forward ); // update 2
normalize( vec, n, forward );
}
} // forwardTrans
/**
Copyright and Use
This software was written and is copyrighted by Ian Kaplan, Bear
Products International, www.bearcave.com, 2001.
iank@bearcave.com
@author Ian Kaplan
*/
public class haar extends liftbase {
/**
Haar predict step
*/
protected void predict( double[] vec, int N, int 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] = vec[j] - predictVal;
}
else if (direction == inverse) {
vec[j] = vec[j] + predictVal;
}
else {
System.out.println("haar::predict: bad direction value");
}
}
}
/**
dj+1, i = oddj+1, i = oddj, i - evenj, i
aj+1, i = evenj, i = (evenj, i + oddj, i)/2
aj+1, i = evenj, i = evenj, i + (oddj, i/2)
*/
protected void update( double[] vec, int N, int 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 {
System.out.println("update: bad direction value");
}
}
}
} // haar
lift/haarpoly.java 100644 1040 1001 15321 7427044775 12733 0 ustar everyone
package lift;
/**
http://www.bearcave.com/misl/misl_tech/wavelets/lifting/index.html
You may use this source code without limitation and without fee as long as you include:
This software was written and is copyrighted by Ian Kaplan, Bear Products International, www.bearcave.com, 2001.
This software is provided "as is", without any warrenty or claim as to its usefulness. Anyone who uses this source code uses it at their own risk. Nor is any support provided by Ian Kaplan and Bear Products International.
Please send any bug fixes or suggested source changes to:
iank@bearcave.com@author Ian Kaplan */ public class haarpoly extends haar { final static int numPts = 4; private polyinterp fourPt; /** haarpoly class constructor */ public haarpoly() { fourPt = new polyinterp(); } /**
Copy four points or N (which ever is less) data points from vec into d These points are the "known" points used in the polynomial interpolation.
@param vec[] the input data set on which the wavelet is calculated @param d[] an array into which N data points, starting at start are copied. @param N the number of polynomial interpolation points @param start the index in vec from which copying starts */ private void fill( double vec[], double d[], int N, int start ) { int n = numPts; if (n > N) n = N; int end = start + n; int j = 0; for (int i = start; i < end; i++) { d[j] = vec[i]; j++; } } // fill /**Predict an odd point from the even points, using 4-point polynomial interpolation.
The four points used in the polynomial interpolation are the even points. We pretend that these four points are located at the x-coordinates 0,1,2,3. The first odd point interpolated will be located between the first and second even point, at 0.5. The next N-3 points are located at 1.5 (in the middle of the four points). The last two points are located at 2.5 and 3.5. For complete documentation see
http://www.bearcave.com/misl/misl_tech/wavelets/lifting/index.html
The difference between the predicted (interpolated) value and the actual odd value replaces the odd value in the forward transform.
As the recursive steps proceed, N will eventually be 4 and then 2. When N = 4, linear interpolation is used. When N = 2, Haar interpolation is used (the prediction for the odd value is that it is equal to the even value).
@param vec the input data on which the forward or inverse transform is calculated. @param N the area of vec over which the transform is calculated @param direction forward or inverse transform */ protected void interp( double[] vec, int N, int direction ) { int half = N >> 1; double d[] = new double[ numPts ]; int k = 42; for (int i = 0; i < half; i++) { double predictVal; if (i == 0) { if (half == 1) { // e.g., N == 2, and we use Haar interpolation predictVal = vec[0]; } else { fill( vec, d, N, 0 ); predictVal = fourPt.interpPoint( 0.5, half, d ); } } else if (i == 1) { predictVal = fourPt.interpPoint( 1.5, half, d ); } else if (i == half-2) { predictVal = fourPt.interpPoint( 2.5, half, d ); } else if (i == half-1) { predictVal = fourPt.interpPoint( 3.5, half, d ); } else { fill( vec, d, N, i-1); predictVal = fourPt.interpPoint( 1.5, half, d ); } int j = i + half; if (direction == forward) { vec[j] = vec[j] - predictVal; } else if (direction == inverse) { vec[j] = vec[j] + predictVal; } else { System.out.println("poly::predict: bad direction value"); } } } // interp /**Haar transform extened with polynomial interpolation forward transform.
This version of the forwardTrans function overrides the function in the liftbase base class. This function introduces an extra polynomial interpolation stage at the end of the transform.
*/ public void forwardTrans( double[] vec ) { final int N = vec.length; for (int n = N; n > 1; n = n >> 1) { split( vec, n ); predict( vec, n, forward ); update( vec, n, forward ); interp( vec, n, forward ); } // for } // forwardTrans /**Haar transform extened with polynomial interpolation inverse transform.
This version of the inverseTrans function overrides the function in the liftbase base class. This function introduces an inverse polynomial interpolation stage at the start of the inverse transform.
*/ public void inverseTrans( double[] vec ) { final int N = vec.length; for (int n = 2; n <= N; n = n << 1) { interp( vec, n, inverse ); update( vec, n, inverse ); predict( vec, n, inverse ); merge( vec, n ); } } // inverseTrans } // haarpoly lift/liftbase.java 100644 1040 1001 14512 7424651473 12702 0 ustar everyone package lift; import util.*; /**class liftbase: base class for simple Lifting Scheme wavelets using split, predict, update or update, predict, merge steps.
Simple lifting scheme wavelets consist of three steps, a split/merge step, predict step and an update step:
The split step divides the elements in an array so that the even elements are in the first half and the odd elements are in the second half.
The merge step is the inverse of the split step. It takes two regions of an array, an odd region and an even region and merges them into a new region where an even element alternates with an odd element.
The predict step calculates the difference between an odd element and its predicted value based on the even elements. The difference between the predicted value and the actual value replaces the odd element.
The predict step operates on the odd elements. The update step operates on the even element, replacing them with a difference between the predict value and the actual odd element. The update step replaces each even element with an average. The result of the update step becomes the input to the next recursive step in the wavelet calculation.
The split and merge methods are shared by all Lifting Scheme wavelet algorithms. This base class provides the transform and inverse transform methods (forwardTrans and inverseTrans). The predict and update methods are abstract and are defined for a particular Lifting Scheme wavelet sub-class.
References:
You may use this source code without limitation and without fee as long as you include:
This software was written and is copyrighted by Ian Kaplan, Bear Products International, www.bearcave.com, 2001.
This software is provided "as is", without any warrenty or claim as to its usefulness. Anyone who uses this source code uses it at their own risk. Nor is any support provided by Ian Kaplan and Bear Products International.
Please send any bug fixes or suggested source changes to:
iank@bearcave.com@author Ian Kaplan */ public abstract class liftbase { /** "enumeration" for forward wavelet transform */ protected final int forward = 1; /** "enumeration" for inverse wavelet transform */ protected final int inverse = 2; /** Split the vec into even and odd elements, where the even elements are in the first half of the vector and the odd elements are in the second half. */ protected void split( double[] vec, int N ) { int start = 1; int end = N - 1; while (start < end) { for (int i = start; i < end; i = i + 2) { double 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. */ protected void merge( double[] 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) { double 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 */ protected abstract void predict( double[] vec, int N, int 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 */ protected abstract void update( double[] vec, int N, int direction ); /**
Simple wavelet Lifting Scheme forward transform
forwardTrans is passed an array of doubles. The array size must be a power of two. 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...).
*/ public void forwardTrans( double[] vec ) { final int N = vec.length; for (int n = N; n > 1; n = n >> 1) { split( vec, n ); predict( vec, n, forward ); update( vec, n, forward ); } } // forwardTrans /**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.
*/ public void inverseTrans( double[] vec ) { final int N = vec.length; for (int n = 2; n <= N; n = n << 1) { update( vec, n, inverse ); predict( vec, n, inverse ); merge( vec, n ); } } // inverseTrans } // liftbase lift/line.java 100644 1040 1001 20410 7425667264 12040 0 ustar everyone package lift; /**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.
You may use this source code without limitation and without fee as long as you include:
This software was written and is copyrighted by Ian Kaplan, Bear Products International, www.bearcave.com, 2001.
This software is provided "as is", without any warrenty or claim as to its usefulness. Anyone who uses this source code uses it at their own risk. Nor is any support provided by Ian Kaplan and Bear Products International.
Please send any bug fixes or suggested source changes to:
iank@bearcave.com@author Ian Kaplan */ public class line extends liftbase { /**
Calculate an extra "even" value for the line wavelet algorithm at the end of the data series. Here we pretend that the last two values in the data series are at the x-axis coordinates 0 and 1, respectively. We then need to calculate the y-axis value at the x-axis coordinate 2. This point lies on a line running through the points at 0 and 1.
Given two points, x1, y1 and x2, y2, where
x1 = 0 x2 = 1
calculate the point on the line at x3, y3, where
x3 = 2
The "two-point equation" for a line given x1, y1 and x2, y2 is
. y2 - y1 (y - y1) = -------- (x - x1) . x2 - x1
Solving for y
. y2 - y1 y = -------- (x - x1) + y1 . x2 - x1
Since x1 = 0 and x2 = 1
. y2 - y1 y = -------- (x - 0) + y1 . 1 - 0
or
y = (y2 - y1)*x + y1
We're calculating the value at x3 = 2, so
y = 2*y2 - 2*y1 + y1
or
y = 2*y2 - y1*/ private double new_y( double y1, double y2) { double y = 2 * y2 - y1; return y; } /**
Predict phase of line Lifting Scheme wavelet
The predict step attempts to "predict" the value of an odd element from the even elements. The difference between the prediction and the actual element is stored as a wavelet coefficient.
The "predict" step takes place after the split step. The split step will move the odd elements (bj) to the second half of the array, leaving the even elements (ai) in the first half
a0, a1, a1, a3, b0, b1, b2, b2,
The predict step of the line wavelet "predicts" that the odd element will be on a line between two even elements.
bj+1,i = bj,i - (aj,i + aj,i+1)/2
Note that when we get to the end of the data series the odd element is the last element in the data series (remember, wavelet algorithms work on data series with 2n elements). Here we "predict" that the odd element will be on a line that runs through the last two even elements. This can be calculated by assuming that the last two even elements are located at x-axis coordinates 0 and 1, respectively. The odd element will be at 2. The new_y() function is called to do this simple calculation.
*/ protected void predict( double[] vec, int N, int 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 predictVal = new_y( vec[i-1], vec[i] ); } if (direction == forward) { vec[j] = vec[j] - predictVal; } else if (direction == inverse) { vec[j] = vec[j] + predictVal; } else { System.out.println("predictline::predict: bad direction value"); } } } // predict /**The predict phase works on the odd elements in the second half of the array. The update phase works on the even elements in the first half of the array. The update phase attempts to preserve the average. After the update phase is completed the average of the even elements should be approximately the same as the average of the input data set from the previous iteration. The result of the update phase becomes the input for the next iteration.
In a Haar wavelet the average that replaces the even element is calculated as the average of the even element and its associated odd element (e.g., its odd neighbor before the split). This is not possible in the line wavelet since the odd element has been replaced by the difference between the odd element and the mid-point of its two even neighbors. As a result, the odd element cannot be recovered.
The value that is added to the even element to preserve the average is calculated by the equation shown below. This equation is given in Wim Sweldens' journal articles and his tutorial (Building Your Own Wavelets at Home) and in Ripples in Mathematics. A somewhat more complete derivation of this equation is provided in Ripples in Mathematics by A. Jensen and A. la Cour-Harbo, Springer, 2001.
The equation used to calculate the average is shown below for a given iteratin i. Note that the predict phase has already completed, so the odd values belong to iteration i+1.
eveni+1,j = eveni,j op (oddi+1,k-1 + oddi+1,k)/4
There is an edge problem here, when i = 0 and k = N/2 (e.g., there is no k-1 element). We assume that the oddi+1,k-1 is the same as oddk. So for the first element this becomes
(2 * oddk)/4
or
oddk/2*/ protected void update( double[] vec, int N, int 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 { System.out.println("update: bad direction value"); } } // for } } // line lift/poly.java 100644 1040 1001 20223 7427047322 12063 0 ustar everyone package lift; /**
Polynomial wavelets
This wavelet transform uses a polynomial interpolation wavelet (e.g., the function used to calculate the differences). A Haar scaling function (the calculation of the average for the even points) is used.
This wavelet transform uses a two stage version of the lifting scheme. In the "classic" two stage Lifting Scheme wavelet the predict stage preceeds the update stage. Also, the algorithm is absolutely symetric, with only the operators (usually addition and subtraction) interchanged.
The problem with the classic Lifting Scheme transform is that it can be difficult to determine how to calculate the smoothing (scaling) function in the update phase once the predict stage has altered the odd values. This version of the wavelet transform calculates the update stage first and then calculates the predict stage from the modified update values. In this case the predict stage uses 4-point polynomial interpolation using even values that result from the update stage.
In this version of the wavelet transform the update stage is no longer perfectly symetric, since the forward and inverse transform equations differ by more than an addition or subtraction operator. However, this version of the transform produces a better result than the Haar transform extended with a polynomial interpolation stage.
This algorithm was suggested to me from my reading of Wim Sweldens' tutorial Building Your Own Wavelets at Home.
http://www.bearcave.com/misl/misl_tech/wavelets/lifting/index.html
You may use this source code without limitation and without fee as long as you include:
This software was written and is copyrighted by Ian Kaplan, Bear Products International, www.bearcave.com, 2001.
This software is provided "as is", without any warrenty or claim as to its usefulness. Anyone who uses this source code uses it at their own risk. Nor is any support provided by Ian Kaplan and Bear Products International.
Please send any bug fixes or suggested source changes to:
iank@bearcave.com@author Ian Kaplan */ public class poly extends liftbase { final static int numPts = 4; private polyinterp fourPt; /** poly class constructor */ public poly() { fourPt = new polyinterp(); } /**
Copy four points or N (which ever is less) data points from vec into d These points are the "known" points used in the polynomial interpolation.
@param vec[] the input data set on which the wavelet is calculated @param d[] an array into which N data points, starting at start are copied. @param N the number of polynomial interpolation points @param start the index in vec from which copying starts */ private void fill( double vec[], double d[], int N, int start ) { int n = numPts; if (n > N) n = N; int end = start + n; int j = 0; for (int i = start; i < end; i++) { d[j] = vec[i]; j++; } } // fill /**The update stage calculates the forward and inverse Haar scaling functions. The forward Haar scaling function is simply the average of the even and odd elements. The inverse function is found by simple algebraic manipulation, solving for the even element given the average and the odd element.
In this version of the wavelet transform the update stage preceeds the predict stage in the forward transform. In the inverse transform the predict stage preceeds the update stage, reversing the calculation on the odd elements.
*/ protected void update( double[] vec, int N, int 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] + vec[j])/2; } else if (direction == inverse) { vec[i] = (2 * vec[i]) - vec[j]; } else { System.out.println("update: bad direction value"); } } } /**Predict an odd point from the even points, using 4-point polynomial interpolation.
The four points used in the polynomial interpolation are the even points. We pretend that these four points are located at the x-coordinates 0,1,2,3. The first odd point interpolated will be located between the first and second even point, at 0.5. The next N-3 points are located at 1.5 (in the middle of the four points). The last two points are located at 2.5 and 3.5. For complete documentation see
http://www.bearcave.com/misl/misl_tech/wavelets/lifting/index.html
The difference between the predicted (interpolated) value and the actual odd value replaces the odd value in the forward transform.
As the recursive steps proceed, N will eventually be 4 and then 2. When N = 4, linear interpolation is used. When N = 2, Haar interpolation is used (the prediction for the odd value is that it is equal to the even value).
@param vec the input data on which the forward or inverse transform is calculated. @param N the area of vec over which the transform is calculated @param direction forward or inverse transform */ protected void predict( double[] vec, int N, int direction ) { int half = N >> 1; double d[] = new double[ numPts ]; int k = 42; for (int i = 0; i < half; i++) { double predictVal; if (i == 0) { if (half == 1) { // e.g., N == 2, and we use Haar interpolation predictVal = vec[0]; } else { fill( vec, d, N, 0 ); predictVal = fourPt.interpPoint( 0.5, half, d ); } } else if (i == 1) { predictVal = fourPt.interpPoint( 1.5, half, d ); } else if (i == half-2) { predictVal = fourPt.interpPoint( 2.5, half, d ); } else if (i == half-1) { predictVal = fourPt.interpPoint( 3.5, half, d ); } else { fill( vec, d, N, i-1); predictVal = fourPt.interpPoint( 1.5, half, d ); } int j = i + half; if (direction == forward) { vec[j] = vec[j] - predictVal; } else if (direction == inverse) { vec[j] = vec[j] + predictVal; } else { System.out.println("poly::predict: bad direction value"); } } } // predict /**Polynomial wavelet lifting scheme transform.
This version of the forwardTrans function overrides the function in the liftbase base class. This function introduces an extra polynomial interpolation stage at the end of the transform.
*/ public void forwardTrans( double[] vec ) { final int N = vec.length; for (int n = N; n > 1; n = n >> 1) { split( vec, n ); update( vec, n, forward ); predict( vec, n, forward ); } // for } // forwardTrans /**Polynomial wavelet lifting Scheme inverse transform.
This version of the inverseTrans function overrides the function in the liftbase base class. This function introduces an inverse polynomial interpolation stage at the start of the inverse transform.
*/ public void inverseTrans( double[] vec ) { final int N = vec.length; for (int n = 2; n <= N; n = n << 1) { predict( vec, n, inverse ); update( vec, n, inverse ); merge( vec, n ); } } // inverseTrans } // poly lift/polyinterp.java 100644 1040 1001 11177 7426421347 13317 0 ustar everyone package lift; class polyinterp { /** number of polynomial interpolation ponts */ private final static int numPts = 4; /** Table for 4-point interpolation coefficients */ private double fourPointTable[][]; /** Table for 2-point interpolation coefficients */ private double twoPointTable[][]; /**The polynomial interpolation algorithm assumes that the known points are located at x-coordinates 0, 1,.. N-1. An interpolated point is calculated at x, using N coefficients. The polynomial coefficients for the point x can be calculated staticly, using the Lagrange method.
@param x the x-coordinate of the interpolated point @param N the number of polynomial points. @param c[] an array for returning the coefficients */ private void lagrange( double x, int N, double c[] ) { double num, denom; for (int i = 0; i < N; i++) { num = 1; denom = 1; for (int k = 0; k < N; k++) { if (i != k) { num = num * (x - k); denom = denom * (i - k); } } // for k c[i] = num / denom; } // for i } // lagrange /**For a given N-point polynomial interpolation, fill the coefficient table, for points 0.5 ... (N-0.5).
*/ private void fillTable( int N, double table[][] ) { double x; double n = N; int i = 0; for (x = 0.5; x < n; x = x + 1.0) { lagrange( x, N, table[i] ); i++; } } // fillTable /**poly constructor
Build the 4-point and 2-point polynomial coefficient tables.
*/ public polyinterp() { // Fill in the 4-point polynomial interplation table // for the points 0.5, 1.5, 2.5, 3.5 fourPointTable = new double[numPts][numPts]; fillTable( numPts, fourPointTable ); // Fill in the 2-point polynomial interpolation table // for 0.5 and 1.5 twoPointTable = new double[2][2]; fillTable( 2, twoPointTable ); } // poly constructor /** Print an N x N table polynomial coefficient table */ private void printTable( double table[][], int N) { System.out.println(N + "-point interpolation table:"); double x = 0.5; for (int i = 0; i < N; i++) { System.out.print(x + ": "); for (int j = 0; j < N; j++) { System.out.print( table[i][j] ); if (j < N-1) System.out.print(", "); } System.out.println(); x = x + 1.0; } } /** Print the 4-point and 2-point polynomial coefficient tables. */ public void printTables() { printTable( fourPointTable, numPts ); printTable( twoPointTable, 2 ); } // printTables /**For the polynomial interpolation point x-coordinate x, return the associated polynomial interpolation coefficients.
@param x the x-coordinate for the interpolated pont @param n the number of polynomial interpolation points @param c[] an array to return the polynomial coefficients */ private void getCoef( double x, int n, double c[]) { double table[][] = null; int j = (int)x; if (j < 0 || j >= n) { System.out.println("poly::getCoef: n = " + n + ", bad x value"); } if (n == numPts) { table = fourPointTable; } else if (n == 2) { table = twoPointTable; c[2] = 0.0; c[3] = 0.0; } else { System.out.println("poly::getCoef: bad value for N"); } if (table != null) { for (int i = 0; i < n; i++) { c[i] = table[j][i]; } } } // getCoef /**Given four points at the x,y coordinates {0,d0}, {1,d1}, {2,d2}, {3,d3} return the y-coordinate value for the polynomial interpolated point at x.
@param x the x-coordinate for the point to be interpolated @param N the number of interpolation points @param d[] an array containing the y-coordinate values for the known points (which are located at x-coordinates 0..N-1). */ public double interpPoint( double x, int N, double d[]) { double c[] = new double[ numPts ]; double point = 0; int n = numPts; if (N < numPts) n = N; getCoef( x, n, c ); if (n == numPts) { point = c[0]*d[0] + c[1]*d[1] + c[2]*d[2] + c[3]*d[3]; } else if (n == 2) { point = c[0]*d[0] + c[1]*d[1]; } return point; } // interpPoint } // polyinterp