`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