`package be.tarsos.dsp.wavelet.lift;/** * <p> * Polynomial wavelets * </p> * <p> * This wavelet transform uses a polynomial interpolation wavelet (e.g., the * function used to calculate the differences). A HaarWavelet scaling function * (the calculation of the average for the even points) is used. * </p> * <p> * 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. * </p> * <p> * 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. * </p> * * <p> * 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 HaarWavelet transform * extended with a polynomial interpolation stage. * </p> * * <p> * This algorithm was suggested to me from my reading of Wim Sweldens' tutorial * <i>Building Your Own Wavelets at Home</i>. * </p> * * </p> * * <pre> * <a href=" http://www.bearcave.com/misl/misl_tech/wavelets/lifting/index.html"> * http://www.bearcave.com/misl/misl_tech/wavelets/lifting/index.html</a> * </pre> * * <h4> * Copyright and Use</h4> * * <p> * You may use this source code without limitation and without fee as long as * you include: * </p> * <blockquote> This software was written and is copyrighted by Ian Kaplan, Bear * Products International, www.bearcave.com, 2001. </blockquote> * <p> * This software is provided "as is", without any 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. * <p> * Please send any bug fixes or suggested source changes to: * * <pre> * iank@bearcave.com * </pre> * * @author Ian Kaplan */public class PolynomialWavelets extends LiftingSchemeBaseWavelet { final static int numPts = 4; private PolynomialInterpolation fourPt; /** * PolynomialWavelets class constructor */ public PolynomialWavelets() { fourPt = new PolynomialInterpolation(); } /** * <p> * Copy four points or <i>N</i> (which ever is less) data points from * <i>vec</i> into <i>d</i> These points are the "known" points used in the * polynomial interpolation. * </p> * * @param vec * the input data set on which the wavelet is calculated * @param d * an array into which <i>N</i> data points, starting at * <i>start</i> are copied. * @param N * the number of polynomial interpolation points * @param start * the index in <i>vec</i> from which copying starts */ private void fill(float vec[], float 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 /** * <p> * The update stage calculates the forward and inverse HaarWavelet scaling * functions. The forward HaarWavelet 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. * </p> * <p> * 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. * </p> */ protected void update(float[] 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"); } } } /** * <p> * Predict an odd point from the even points, using 4-point polynomial * interpolation. * </p> * <p> * 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 * </p> * * <pre> * <a href=" http://www.bearcave.com/misl/misl_tech/wavelets/lifting/index.html"> * http://www.bearcave.com/misl/misl_tech/wavelets/lifting/index.html</a> * </pre> * * <p> * The difference between the predicted (interpolated) value and the actual * odd value replaces the odd value in the forward transform. * </p> * * <p> * As the recursive steps proceed, N will eventually be 4 and then 2. When N * = 4, linear interpolation is used. When N = 2, HaarWavelet interpolation * is used (the prediction for the odd value is that it is equal to the even * value). * </p> * * @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(float[] vec, int N, int direction) { int half = N >> 1; float d[] = new float[numPts]; // int k = 42; for (int i = 0; i < half; i++) { float predictVal; if (i == 0) { if (half == 1) { // e.g., N == 2, and we use HaarWavelet interpolation predictVal = vec; } else { fill(vec, d, N, 0); predictVal = fourPt.interpPoint(0.5f, half, d); } } else if (i == 1) { predictVal = fourPt.interpPoint(1.5f, half, d); } else if (i == half - 2) { predictVal = fourPt.interpPoint(2.5f, half, d); } else if (i == half - 1) { predictVal = fourPt.interpPoint(3.5f, half, d); } else { fill(vec, d, N, i - 1); predictVal = fourPt.interpPoint(1.5f, 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("PolynomialWavelets::predict: bad direction value"); } } } // predict /** * <p> * Polynomial wavelet lifting scheme transform. * </p> * <p> * This version of the forwardTrans function overrides the function in the * LiftingSchemeBaseWavelet base class. This function introduces an extra * polynomial interpolation stage at the end of the transform. * </p> */ public void forwardTrans(float[] 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 /** * <p> * Polynomial wavelet lifting Scheme inverse transform. * </p> * <p> * This version of the inverseTrans function overrides the function in the * LiftingSchemeBaseWavelet base class. This function introduces an inverse * polynomial interpolation stage at the start of the inverse transform. * </p> */ public void inverseTrans(float[] 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} // PolynomialWavelets`