This is my programmatic note 8 on "Ripples in Mathematics"
by A. Jensen & A. la Cour-Harbo. This note is on Chapters 3 and 4 where CDF(2, 2) and CDF(4, 4). As I mentioned in my programmatic note 7 on Ripples, I ran into problems trying to reconstruct the graphs in Fig. 4.10 and 4.11 in Ch.
4 with CDF(2,2). My graphs were not nearly as smooth as the graphs in the book. When I dug deeper I discovered that the forward and inverse normalization coefficients for CDF(2, 2) varied from publication to publication. So I decided to switch to CDF(4, 4), because the normalization coefficients seem to agree across publications.
The multiresolution analysis described by Jensen and la Cour-Harbo requires computing forward and inverse DWT for a given number of iterations or scales. This blog entry is on computing CDF(4, 4) for a given number of iterations. The constants H0, H1, H2, H3, G0, G1, G2, and G2 are defined in my note 7. My Java source is here. In the next two programmatic notes on Ripples, I hope to address the computation of inverse CDF(4, 4) for all scales and for a given number of scales.
The multiresolution analysis described by Jensen and la Cour-Harbo requires computing forward and inverse DWT for a given number of iterations or scales. This blog entry is on computing CDF(4, 4) for a given number of iterations. The constants H0, H1, H2, H3, G0, G1, G2, and G2 are defined in my note 7. My Java source is here. In the next two programmatic notes on Ripples, I hope to address the computation of inverse CDF(4, 4) for all scales and for a given number of scales.
I started by defining constants for H0, H1, H2, H3 and G0, G1, G2, and G3. The second parameter is the number of iterations we want CDF(4, 4) to run on the signal. For example, if the signal's length is 4, we can run 1 iteration, because CDF(4, 4) stops when the signal length is 4. When the signal's length is 16, we can run 2 iterations, etc. The third parameter dbg_flag variable that can be set to false when output is not needed. Utils.isPowerOf2(N) is from my Utils class. It has a bunch of util methods, primarily numerical.
public static void orderedDWTForNumIters(double[] signal, int num_iters, boolean dbg_flag) {
final int N = signal.length;
if ( N < 4 || !Utils.isPowerOf2(N) ) {
System.out.println("No DWT will be done: signal's length is < 4 or not a power of 2");
return;
}
int i, j, mid;
double[] D4 = null;
if ( dbg_flag ) { System.out.print("=>INPUT: "); Utils.displaySample(signal); }
int numScalesToDo = Utils.powVal(N)-1;
int currScale = 0;
int signal_length = N;
while ( signal_length >= 4 ) {
mid = signal_length >> 1; // n / 2;
if ( dbg_flag ) System.out.println("MID = " + mid);
if ( dbg_flag ) System.out.println("signal_length = " + signal_length);
D4 = new double[signal_length]; // temporary array that saves the scalers and wavelets
for(i = 0, j = 0; j < signal_length-3; i += 1, j += 2) {
if ( dbg_flag ) {
final String cursig = "s^{" + (currScale+1) + "}_{" + (numScalesToDo-1) + "}";
final String prvsig = "s^{" + currScale + "}_{" + numScalesToDo + "}";
System.out.println("FWD SCL: " + cursig + "[" + i + "]=" + "H0*" + prvsig + "[" + j + "]+H1*" + prvsig + "[" + (j+1) + "]+" +
"H2*" + prvsig + "[" + (j+2) + "]+" + "H3*" + prvsig + "[" + (j+3) + "]; " );
System.out.println("FWD WVL: " + cursig + "[" + (mid+i) + "]=" + "G0*" + prvsig + "[" + j + "]+" + "G1*" + prvsig + "[" + (j+1) + "]+" +
"G2*" + prvsig + "[" + (j+2) + "]+" + "G3*" + prvsig + "[" + (j+3) + "]" );
}
// cdf44[i] is a scaled sample
D4[i] = H0*signal[j] + H1*signal[j+1] + H2*signal[j+2] + H3*signal[j+3];
// cdf44[mid+i] is the corresponding wavelet for d4[i]
D4[mid+i] = G0*signal[j] + G1*signal[j+1] + G2*signal[j+2] + G3*signal[j+3];
}
// cdf44[i] is a scaled sample with a mirror wrap-up
D4[i] = H0*signal[signal_length-2] + H1*signal[signal_length-1] + H2*signal[0] + H3*signal[1];
// cdf44[mid+i] is the corresponding wavelet for d4[i]
D4[mid+i] = G0*signal[signal_length-2] + G1*signal[signal_length-1] + G2*signal[0] + G3*signal[1];
if ( dbg_flag ) {
final String cursig = "s^{" + currScale + "}_{" + numScalesToDo + "}";
final String prvsig = "s^{" + (currScale-1) + "}_{" + (numScalesToDo+1) + "}";
System.out.println("FWD SCL: " + cursig + "[" + i + "]=" + "H0*" + prvsig + "[" + (signal_length-2) + "]+H1*" +
public static void orderedDWTForNumIters(double[] signal, int num_iters, boolean dbg_flag) {
final int N = signal.length;
if ( N < 4 || !Utils.isPowerOf2(N) ) {
System.out.println("No DWT will be done: signal's length is < 4 or not a power of 2");
return;
}
int i, j, mid;
double[] D4 = null;
if ( dbg_flag ) { System.out.print("=>INPUT: "); Utils.displaySample(signal); }
int numScalesToDo = Utils.powVal(N)-1;
int currScale = 0;
int signal_length = N;
while ( signal_length >= 4 ) {
mid = signal_length >> 1; // n / 2;
if ( dbg_flag ) System.out.println("MID = " + mid);
if ( dbg_flag ) System.out.println("signal_length = " + signal_length);
D4 = new double[signal_length]; // temporary array that saves the scalers and wavelets
for(i = 0, j = 0; j < signal_length-3; i += 1, j += 2) {
if ( dbg_flag ) {
final String cursig = "s^{" + (currScale+1) + "}_{" + (numScalesToDo-1) + "}";
final String prvsig = "s^{" + currScale + "}_{" + numScalesToDo + "}";
System.out.println("FWD SCL: " + cursig + "[" + i + "]=" + "H0*" + prvsig + "[" + j + "]+H1*" + prvsig + "[" + (j+1) + "]+" +
"H2*" + prvsig + "[" + (j+2) + "]+" + "H3*" + prvsig + "[" + (j+3) + "]; " );
System.out.println("FWD WVL: " + cursig + "[" + (mid+i) + "]=" + "G0*" + prvsig + "[" + j + "]+" + "G1*" + prvsig + "[" + (j+1) + "]+" +
"G2*" + prvsig + "[" + (j+2) + "]+" + "G3*" + prvsig + "[" + (j+3) + "]" );
}
// cdf44[i] is a scaled sample
D4[i] = H0*signal[j] + H1*signal[j+1] + H2*signal[j+2] + H3*signal[j+3];
// cdf44[mid+i] is the corresponding wavelet for d4[i]
D4[mid+i] = G0*signal[j] + G1*signal[j+1] + G2*signal[j+2] + G3*signal[j+3];
}
// cdf44[i] is a scaled sample with a mirror wrap-up
D4[i] = H0*signal[signal_length-2] + H1*signal[signal_length-1] + H2*signal[0] + H3*signal[1];
// cdf44[mid+i] is the corresponding wavelet for d4[i]
D4[mid+i] = G0*signal[signal_length-2] + G1*signal[signal_length-1] + G2*signal[0] + G3*signal[1];
if ( dbg_flag ) {
final String cursig = "s^{" + currScale + "}_{" + numScalesToDo + "}";
final String prvsig = "s^{" + (currScale-1) + "}_{" + (numScalesToDo+1) + "}";
System.out.println("FWD SCL: " + cursig + "[" + i + "]=" + "H0*" + prvsig + "[" + (signal_length-2) + "]+H1*" +
prvsig + "[" + (signal_length-1) + "]+" +
"H2*" + prvsig + "[" + 0 + "]+" + "H3*" + prvsig + "[" + 1 + "]; " );
System.out.println("FWD WVL: " + cursig + "[" + (mid+i) + "]=" + "G0*" + prvsig + "[" + (signal_length-2) + "]+"
"H2*" + prvsig + "[" + 0 + "]+" + "H3*" + prvsig + "[" + 1 + "]; " );
System.out.println("FWD WVL: " + cursig + "[" + (mid+i) + "]=" + "G0*" + prvsig + "[" + (signal_length-2) + "]+"
+ "G1*" + prvsig + "[" + (signal_length-1) + "]+" +
"G2*" + prvsig + "[" + 0 + "]+" + "G3*" + prvsig + "[" + 1 + "]" );
}
System.arraycopy(D4, 0, signal, 0, D4.length);
D4 = null;
signal_length >>= 1; // signal_length gets halved at each iteration/scale
currScale += 1;
numScalesToDo -= 1;
if ( currScale >= num_iters ) return;
}
}
"G2*" + prvsig + "[" + 0 + "]+" + "G3*" + prvsig + "[" + 1 + "]" );
}
System.arraycopy(D4, 0, signal, 0, D4.length);
D4 = null;
signal_length >>= 1; // signal_length gets halved at each iteration/scale
currScale += 1;
numScalesToDo -= 1;
if ( currScale >= num_iters ) return;
}
}
Below are some sample signals, a test method, and a couple tests in the main with the dbg_flag set to false.
static double[] a01a = {1, 2, 3, 4};
static double[] a01b = {4, 3, 2, 1};
static double[] a02a = {1, 2, 3, 4, 5, 6, 7, 8};
static double[] a02b = {8, 7, 6, 5, 4, 3, 2, 1};
static double[] a03a = {1, 1, 1, 1};
static double[] a03b = {2, 2, 2, 2};
static double[] a03c = {3, 3, 3, 3};
static double[] a03d = {4, 4, 4, 4};
static double[] a04a = {1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16};
static double[] a04b = {16, 15, 14, 13, 12, 11, 10, 9, 8, 7, 6, 5, 4, 3, 2, 1};
public static void main(String[] args) {
test_fwd_inv_cdf44(a02a, false);
test_fwd_inv_cdf44(a02b, false);
}
Here is the output of the main.
Input: 1.0 2.0 3.0 4.0 5.0 6.0 7.0 8.0
FWD CDF(4,4): 5.901923788646682 12.098076211353314 0.3660254037844384 -3.8301270189221923 -4.440892098500626E-16 -8.881784197001252E-16 -8.881784197001252E-16 -2.82842712474619
INV CDF(4,4): 1.0 1.9999999999999998 2.999999999999998 3.999999999999997 4.9999999999999964 5.999999999999997 6.999999999999997 7.999999999999998
Input: 8.0 7.0 6.0 5.0 4.0 3.0 2.0 1.0
FWD CDF(4,4): 12.098076211353314 5.901923788646685 -0.3660254037844395 3.8301270189221936 -8.881784197001252E-16 -4.440892098500626E-16 -2.220446049250313E-16 2.828427124746189
INV CDF(4,4): 7.9999999999999964 6.999999999999998 5.999999999999997 4.999999999999997 3.9999999999999987 2.9999999999999996 2.0 0.9999999999999997
static double[] a01b = {4, 3, 2, 1};
static double[] a02a = {1, 2, 3, 4, 5, 6, 7, 8};
static double[] a02b = {8, 7, 6, 5, 4, 3, 2, 1};
static double[] a03a = {1, 1, 1, 1};
static double[] a03b = {2, 2, 2, 2};
static double[] a03c = {3, 3, 3, 3};
static double[] a03d = {4, 4, 4, 4};
static double[] a04a = {1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16};
static double[] a04b = {16, 15, 14, 13, 12, 11, 10, 9, 8, 7, 6, 5, 4, 3, 2, 1};
public static void main(String[] args) {
test_fwd_inv_cdf44(a02a, false);
test_fwd_inv_cdf44(a02b, false);
}
Here is the output of the main.
Input: 1.0 2.0 3.0 4.0 5.0 6.0 7.0 8.0
FWD CDF(4,4): 5.901923788646682 12.098076211353314 0.3660254037844384 -3.8301270189221923 -4.440892098500626E-16 -8.881784197001252E-16 -8.881784197001252E-16 -2.82842712474619
INV CDF(4,4): 1.0 1.9999999999999998 2.999999999999998 3.999999999999997 4.9999999999999964 5.999999999999997 6.999999999999997 7.999999999999998
Input: 8.0 7.0 6.0 5.0 4.0 3.0 2.0 1.0
FWD CDF(4,4): 12.098076211353314 5.901923788646685 -0.3660254037844395 3.8301270189221936 -8.881784197001252E-16 -4.440892098500626E-16 -2.220446049250313E-16 2.828427124746189
INV CDF(4,4): 7.9999999999999964 6.999999999999998 5.999999999999997 4.999999999999997 3.9999999999999987 2.9999999999999996 2.0 0.9999999999999997
