plectrum

PolynomialWavelets.java (8436B)

---

 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[0];
 				} 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