lifttest.java100644 1040 1001 4075 7427050360 11763 0ustar 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.cpp100644 1040 1001 1226 7426143205 11617 0ustar everyone #include #include using namespace std; typedef vector DoubleVec; #include "C++\lifting.h" 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 }; void main() { debug dbg; liftingScheme poly; int N = sizeof( data ) / sizeof( double ); DoubleVec ts; int i; for (i = 0; i < N; i++) { ts.push_back( data[i] ); } printf("data: "); dbg.pr_vector( ts ); printf("\n"); poly.transform( ts ); dbg.pr_ordered( ts ); printf("\n"); poly.invTransform( ts ); dbg.pr_vector( ts ); printf("\n"); } C++/ 40755 1040 1001 0 7427050620 7466 5ustar everyoneC++/interp.h100644 1040 1001 14135 7366611142 11265 0ustar everyone /** Support for four point polynomial interpolation, using the Lagrange formula. */ class interpolation { private: double coef_minus2[4]; double coef_minus1[4]; /** lagrange coefficients for x = -0.5 */ double coef_minus05[4]; /** lagrange coefficients for x = 0.5 */ double coef_05[4]; /** lagrange coefficients for x = 1.5 */ double coef_15[4]; /** lagrange coefficients for x = 2.5 */ double coef_25[4]; /** lagrange coefficients for x = 3.5 */ double coef_35[4]; double coef_40[4]; double coef_50[4]; /** Evaluate the Lagrange formula at x producing a set of coefficients for the equation 
      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.h100644 1040 1001 3055 7363726044 11424 0ustar 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.h100644 1040 1001 43246 7426422112 11417 0ustar everyone #include /** \file This file contains code to implement Lifting Scheme wavelets. Lifting Scheme wavelets are described in Wim Sweldens' tutorial paper Building Your Own Wavelets at Home which is available on the Web. Lifting Scheme wavelets are a conceptual extension of Haar wavelets. One of the disadvantages of Haar wavelets is the high frequency (largest) coefficient spectrum can miss detail (even to odd transitions, for example). Lifting Scheme wavelets properly represent change in all coefficient spectrum. This makes lifting scheme wavelets a better choice for some algorithms which do filtering based on the absolute standard deviation calculated on the high frequency coefficients. */ /** This class contains functions that are useful in debugging wavelet code. */ class debug { public: /** Print the result of a wavelet transform so that each coefficient spectrum is printed in its own block, showing the structure of the result. */ void pr_ordered( const DoubleVec& coef) { const char *format = "%7.4f "; printf("{"); int len = coef.size(); int cnt = 0; int num_in_freq = 0; for (int i = 0; i < len; i++) { printf(format, coef[i]); cnt++; if (num_in_freq == 0) { printf("\n"); cnt = 0; num_in_freq = 1; } else if (cnt == num_in_freq) { printf("\n"); cnt = 0; num_in_freq = num_in_freq * 2; } } printf("}\n"); } // pr_ordered /** Print wavelet transform intermediate results. The wavelet transform calculates a set of averages and coefficients (differences). These are calculated recursively. For example, for a time series with size 32, the first iteration of the recursive wavelet transform will calculate a set of 16 averages and 16 coefficients (differences). The next set of differences and avarages will be calculated at the start of the array (0..15) and will consist of 8 averages and 8 differences. This function is passed a value for N which is the size of the average and difference spectrum. It prints each of these. */ void pr_intermediate( DoubleVec& ts, int N) { const char* fmt = "%7.4f "; int i; int half = N >> 1; printf(" ave = "); for (i = 0; i < half; i++) { printf(fmt, ts[i] ); } printf("\n"); printf("coef = "); for (i = half; i < N; i++) { printf(fmt, ts[i] ); } printf("\n"); } // pr_intermediate /** Print the first N elements of ts. If N is less than the length of ts the print a dividing mark (e.g., "|") and print the rest of the array. This function is useful for displaying the flow of the wavelet computation. */ void pr_vector( DoubleVec& ts) { const int N = ts.size(); const char* fmt = "%7.4f "; int i; for (i = 0; i < N; i++) { printf(fmt, ts[i] ); } int len = ts.size(); if (len != N) { printf(" | "); for (i = N; i < len; i++) { printf(fmt, ts[i] ); } } printf("\n"); } // pr_vector }; // debug /** Support for four point polynomial interpolation, using the Lagrange formula. */ class interpolation { private: double fourPointTable[4][4]; double twoPointTable[4][4]; /**

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 */ 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).

*/ void fillTable( const int N, double table[4][4] ) { 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 /** Print an N x N table polynomial coefficient table */ void printTable( double table[4][4], int N) { printf("%d-point interpolation table:\n", N); double x = 0.5; for (int i = 0; i < N; i++) { printf("%4.2f: ", x); for (int j = 0; j < N; j++) { printf("%6.4f", table[i][j] ); if (j < N-1) printf(", "); } printf("\n"); x = x + 1.0; } } /** Print the 4-point and 2-point polynomial coefficient tables. */ void printTables() { printTable( fourPointTable, 4 ); printTable( twoPointTable, 2 ); printf("\n"); } // 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 */ void getCoef( double x, int n, double c[]) { int j = (int)x; if (j < 0 || j >= n) { printf("poly::getCoef: n = %d, bad x value = %f\n", n, x); } if (n == 2) { c[2] = 0.0; c[3] = 0.0; } else if (n != 4) { printf("poly::getCoef: bad value for N\n"); return; } for (int i = 0; i < n; i++) { if (n == 4) { c[i] = fourPointTable[j][i]; } else { // n == 2 c[i] = twoPointTable[j][i]; } } } // getCoef 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() { // Fill in the 4-point polynomial interplation table // for the points 0.5, 1.5, 2.5, 3.5 fillTable( 4, fourPointTable ); // Fill in the 2-point polynomial interpolation table // for 0.5 and 1.5 fillTable( 2, twoPointTable ); } // interpolation class constructor /**

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). */ double interpPoint( double x, int N, double d[]) { double c[4]; double point = 0; int n = 4; if (N < 4) n = N; getCoef( x, n, c ); if (n == 4) { 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 }; // class interpolation /** This class implements the wavelet Lifting Scheme with polynomial interpolation. This version of wavelets is based on Wim Sweldens' tutorial paper Building Your Own Wavelets at Home. This tutorial was presented at SIGGraph. The tutorial is available on the Web. One of the attractive features of Lifting Scheme wavelets is that the transform and the inverse transform are exact mirrors of each other. To arrive at the inverse transform the order of the operations is reversed and subtraction and addition operations are exchanged. This allows a generic Lifting Scheme base class to be constructed, where the ordering is defined by the sub-classes. */ class liftingScheme { private: interpolation fourPtInterp; /** Enumerator for add and subtract operators */ typedef enum { bad_dir, forward, inverse } direction; /**

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 */ void fill( DoubleVec& vec, double d[], int N, int start ) { int n = 4; 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 */ void interp( DoubleVec& vec, const int N, const direction dir ) { int half = N >> 1; double d[4]; 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 = fourPtInterp.interpPoint( 0.5, half, d ); } } else if (i == 1) { predictVal = fourPtInterp.interpPoint( 1.5, half, d ); } else if (i == half-2) { predictVal = fourPtInterp.interpPoint( 2.5, half, d ); } else if (i == half-1) { predictVal = fourPtInterp.interpPoint( 3.5, half, d ); } else { fill( vec, d, N, i-1); predictVal = fourPtInterp.interpPoint( 1.5, half, d ); } int j = i + half; if (dir == forward) { vec[j] = vec[j] - predictVal; } else if (dir == inverse) { vec[j] = vec[j] + predictVal; } else { printf("interp: bad direction value\n"); break; } } // for } // interp /** For the first N elements of ts, move the odd elements to the second half of the vector, leaving the even elements in the first half. In terms of the wavelet transform, this divides the vector into an average region and a coefficient (difference) region. For example, if we label the even elements with a and the odd elements with b the input vector will contain 
  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 5ustar everyonelift/daubechies.java100644 1040 1001 6761 7424663416 13174 0ustar 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 /**

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) { normalize( vec, n, inverse ); update( vec, n, inverse ); predict( vec, n, inverse ); updateOne(vec, n, inverse ); merge( vec, n ); } } // inverseTrans } lift/haar.java100644 1040 1001 7667 7425667444 12027 0ustar everyone package lift; /**

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, ...)

Copyright and Use

You may use this source code without limitation and without fee as long as you include:

This software was written and is copyrighted by Ian Kaplan, Bear Products International, www.bearcave.com, 2001.

This software is provided "as is", without any 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 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"); } } } /**

Update step of the Haar wavelet transform.

The wavelet transform calculates a set of detail or difference coefficients in the predict step. These are stored in the upper half of the array. The update step calculates an average from the even-odd element pairs. The averages will replace the even elements in the lower half of the array.

The Haar wavelet calculation used in the Lifting Scheme is 

       dj+1, i = oddj+1, i = oddj, i - evenj, i
       aj+1, i = evenj, i = (evenj, i + oddj, i)/2

Note that the Lifting Scheme uses an in-place algorithm. The odd elements have been replaced by the detail coefficients in the predict step. With a little algebra we can substitute the coefficient calculation into the average calculation, which gives us 

       aj+1, i = evenj, i = evenj, i + (oddj, i/2)
*/ 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.java100644 1040 1001 15321 7427044775 12733 0ustar everyone package lift; /**

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. 

  
  http://www.bearcave.com/misl/misl_tech/wavelets/lifting/index.html

Copyright and Use

You may use this source code without limitation and without fee as long as you include:

This software was written and is copyrighted by Ian Kaplan, Bear Products International, www.bearcave.com, 2001.

This software is provided "as is", without any 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.java100644 1040 1001 14512 7424651473 12702 0ustar 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:

  • The Wavelet Lifting Scheme by Ian Kaplan, www.bearcave.com. This is the parent web page for this Java source code.
  • Ripples in Mathematics: the Discrete Wavelet Transform by Arne Jense and Anders la Cour-Harbo, Springer, 2001
  • Building Your Own Wavelets at Home in Wavelets in Computer Graphics

Copyright and Use

You may use this source code without limitation and without fee as long as you include:

This software was written and is copyrighted by Ian Kaplan, Bear Products International, www.bearcave.com, 2001.

This software is provided "as is", without any 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.java100644 1040 1001 20410 7425667264 12040 0ustar 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.

Copyright and Use

You may use this source code without limitation and without fee as long as you include:

This software was written and is copyrighted by Ian Kaplan, Bear Products International, www.bearcave.com, 2001.

This software is provided "as is", without any 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.java100644 1040 1001 20223 7427047322 12063 0ustar 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

Copyright and Use

You may use this source code without limitation and without fee as long as you include:

This software was written and is copyrighted by Ian Kaplan, Bear Products International, www.bearcave.com, 2001.

This software is provided "as is", without any 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.java100644 1040 1001 11177 7426421347 13317 0ustar 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