The book Ripples in Mathematics by A. Jensen and A. la Cour-Harbo, Springer Verlag, 2001, includes a C version of the Lifting Scheme Daubechies D4 forward transform. Perhaps to show the correspondence between this code and the equations that describe the lifting scheme version of the Daubechies D4 transform, the code uses temporaries (vectors d and s). This code is reproduced here without these temporaries. Also an inverse transform step is included.
The lifting scheme divides wavelet transforms into a set of steps. In other cases I've implemented each step as a function. In this case, following Ripples, all of the steps for the forward and inverse transform are included in a single function.
Rather than attempt to reproduce the background of the lifting scheme version of the Daubechies D4 transform here, I will refer the reader to Ripples.
author: Ian Kaplan, January 2002
*/ class daubbook { final static double sqrt3 = Math.sqrt( 3 ); final static double sqrt2 = Math.sqrt( 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; } } /** Forward step of the Daubechies D4 transform */ protected void forwardStep( double[] S, int N ) { final int half = N/2; // update 1 for (int n = 0; n < half; n++) S[n] = S[n] + sqrt3 * S[half+n]; // predict S[half] = S[half] - (sqrt3/4.0)*S[0] - (((sqrt3-2)/4.0)*S[half-1]); for (int n = 1; n < half; n++) S[half+n] = S[half+n] - (sqrt3/4.0)*S[n] - (((sqrt3-2)/4.0)*S[n-1]); // update 2 for (int n = 0; n < half-1; n++) S[n] = S[n] - S[half+n+1]; S[half-1] = S[half-1] - S[half]; // normalize for (int n = 0; n < half; n++) { S[n] = ((sqrt3-1.0)/sqrt2) * S[n]; S[n+half] = ((sqrt3+1.0)/sqrt2) * S[n+half]; } } // forwardStep /** Inverse step of the Daubechies D4 transform */ protected void inverseStep( double[] S, int N ) { final int half = N/2; for (int n = 0; n < half; n++) { S[n] = ((sqrt3+1.0)/sqrt2) * S[n]; S[n+half] = ((sqrt3-1.0)/sqrt2) * S[n+half]; } for (int n = 0; n < half-1; n++) S[n] = S[n] + S[half+n+1]; S[half-1] = S[half-1] + S[half]; S[half] = S[half] + (sqrt3/4.0)*S[0] + (((sqrt3-2)/4.0)*S[half-1]); for (int n = 1; n < half; n++) S[half+n] = S[half+n] + (sqrt3/4.0)*S[n] + (((sqrt3-2)/4.0)*S[n-1]); for (int n = 0; n < half; n++) S[n] = S[n] - sqrt3 * S[half+n]; } // inverseStep /** Forward Daubechies transform, implemented with a single monolithic forward step. */ public void forwardTrans( double[] vec ) { final int N = vec.length; for (int n = N; n > 1; n = n >> 1) { split(vec, n ); forwardStep( vec, n ); } } // forwardTrans /** Inverse Daubechies transform, implemented with a single monolithic inverse step. */ public void inverseTrans( double[] vec ) { final int N = vec.length; for (int n = 2; n <= N; n = n << 1) { inverseStep( vec, n ); merge(vec, n ); } } // inverseTrans } // daubbook lift/liftbase.java 100644 1040 1001 14512 7424651473 12702 0 ustar 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:
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/daubechies.java 100644 1040 1001 6761 7424663416 13174 0 ustar 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 }