Daubechies4Wavelet.java

package be.tarsos.dsp.wavelet.lift;

/**
* 
* @author Ian Kaplan
*/
public class Daubechies4Wavelet extends LiftingSchemeBaseWavelet {
	final static float sqrt3 = (float) Math.sqrt(3);
	final static float sqrt2 = (float) Math.sqrt(2);

	protected void normalize(float[] S, int N, int direction) {
		int half = N >> 1;

		for (int n = 0; n < half; n++) {
			if (direction == forward) {
				S[n] = ((sqrt3 - 1.0f) / sqrt2) * S[n];
				S[n + half] = ((sqrt3 + 1.0f) / sqrt2) * S[n + half];
			} else if (direction == inverse) {
				S[n] = ((sqrt3 + 1.0f) / sqrt2) * S[n];
				S[n + half] = ((sqrt3 - 1.0f) / sqrt2) * S[n + half];
			} else {
				System.out
						.println("Daubechies4Wavelet::normalize: bad direction value");
				break;
			}
		}
	} // normalize

	protected void predict(float[] S, int N, int direction) {
		int half = N >> 1;

		if (direction == forward) {
			S[half] = S[half] - (sqrt3 / 4.0f) * S[0]
					- (((sqrt3 - 2) / 4.0f) * S[half - 1]);
		} else if (direction == inverse) {
			S[half] = S[half] + (sqrt3 / 4.0f) * S[0]
					+ (((sqrt3 - 2) / 4.0f) * S[half - 1]);
		} else {
			System.out
					.println("Daubechies4Wavelet::predict: bad direction value");
		}

		// predict, forward

		for (int n = 1; n < half; n++) {
			if (direction == forward) {
				S[half + n] = S[half + n] - (sqrt3 / 4.0f) * S[n]
						- (((sqrt3 - 2) / 4.0f) * S[n - 1]);
			} else if (direction == inverse) {
				S[half + n] = S[half + n] + (sqrt3 / 4.0f) * S[n]
						+ (((sqrt3 - 2) / 4.0f) * S[n - 1]);
			} else {
				break;
			}
		}

	} // predict

	protected void updateOne(float[] S, int N, int direction) {
		int half = N >> 1;

		for (int n = 0; n < half; n++) {
			float 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("Daubechies4Wavelet::updateOne: bad direction value");
				break;
			}
		}
	} // updateOne

	protected void update(float[] 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("Daubechies4Wavelet::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(float[] 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

	/**
	 * <p>
	 * Default two step Lifting Scheme inverse wavelet transform
	 * </p>
	 * 
	 * <p>
	 * 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.
	 * </p>
	 */
	public void inverseTrans(float[] 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
}