Computing Forward Ordered CDF(4, 4): Programmatic Note 7 on Ripples in Mathematics

Vladimir Kulyukin

---

Problem
  

This is my programmatic note 7 on "Ripples in Mathematics" by A. Jensen & A. la Cour-Harbo. One of the problems I had when I was working on CDF(2, 2) was the inverse transform. My values, when I was inversing the transformed signals, were off by 0.25 or some small amounts. I could not reconstruct the graphs in Fig. 4.10 and 4.11 in Ch. 4 in "Ripples in Mathematics." In general, my CDF implementation did not go nearly as smoothly as the Haar. I have documented my conceptual struggles and formalism modifications in my previous posts on CDF(2, 2) on this blog which you can find by searching it on "wavelets." 

I started looking deeper into the coefficients for computing smoothed values and wavelets, i.e., H0, H1, G0, G1 for CDF(2, 2) and H0, H1, H2, H3, G0, G1, G2, G3, G4 for forward CDF(4, 4).  What I found out was that these coefficients vary from publication to publication. Perhaps, I am still off on some of the finer mathematical points of Daubechies Wavelets and simply do not understand a formal abstraction (if one exists) where they are all the same. For those of you who are interested in digging deeper into these coefficients, here are the three references I looked at to better understand CDF(2, 2) and CDF(4, 4). 

1) Y. Nievergelt. "Wavelets Made Easy";
2) G. Uytterhoeven, D. Roose, and A. Buttheel. "Wavelet Transforms Using the Lifting Scheme."; 
3) G. Uytterhoeven, F. Van Wulpen, M. Jansen, D. Roose, and A. Bultheel. "WAILI: Wavelets with Integer Lifting."

These are great references and I keep going back to them and re-reading them, but they are very formal and lack programmatic details, especially for software developers working with imperative programming languages such as Java and C++. As luck would have it, I then found Ian Kaplan's wonderful site at www.bearcave.com. Ian has a whole bunch of great documents and C++ and Java programs on wavelets. After I read his explanation of CDF(4, 4), a lot of things became much clearer to me programmatically. Then  I got back to implementing forward CDF(4, 4) and used Ian's coefficients. He credits Jensen and la Cour-Harbo in his implementation of CDF(4, 4). However, I could not find where in the book those coefficients are defined. I will keep looking for them though. Below is my implementation of forward CDF(4, 4). I plan to post my inverse implementation in a future blog entry. My Java source is here.

---

Forward CDF(4, 4) in Java
  

I started by defining constants for H0, H1, H2, H3 and G0, G1, G2, and G3.

public class CDF44 {
    static final double SQRT_OF_3 = Math.sqrt(3);
    static final double SQRT_OF_2 = Math.sqrt(2);
    static final double FOUR_SQRT_OF_2 = 4*SQRT_OF_2;
   
    // CDF(4,4) Forward signal coefficients
    static final double H0 = (1 + SQRT_OF_3)/FOUR_SQRT_OF_2;
    static final double H1 = (3 + SQRT_OF_3)/FOUR_SQRT_OF_2;
    static final double H2 = (3 - SQRT_OF_3)/FOUR_SQRT_OF_2;
    static final double H3 = (1 - SQRT_OF_3)/FOUR_SQRT_OF_2;
   
    // CDF(4, 4) Forward wavelet coefficients
    static final double G0 = H3;
    static final double G1 = -H2;
    static final double G2 = H1;
    static final double G3 = -H0;

}

Here is an implementation of forward CDF(4, 4) as a static method of CDF44. There is a dbg_flag variable that can be set to false when output is not needed. I generally prefer while-loops to for-loops. There are, and have always been, more intuitively obvious to me. Utils.isPowerOf2(N) is from my Utils class. It has a bunch of util methods, primarily numerical.

    public static void orderedDWT(double[] signal, boolean dbg_flag) {
        final int N = signal.length;
        if ( !Utils.isPowerOf2(N) ) return;
        if ( N < 4 ) return;
        int i, j, mid;
        double[] D4 = null;

        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.print("SCL:  " + cursig + "[" + i + "]=" + "H0*" + prvsig + "[" + j + "]+H1*" + prvsig + "[" + (j+1) + "]+" +
                        "H2*" + prvsig + "[" + (j+2) + "]+" + "H3*" + prvsig + "[" + (j+3) + "]; " );
                    System.out.println("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];
            }

            currScale     += 1;
            numScalesToDo -= 1;
          
            // 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.print("SCL:  " + cursig + "[" + i + "]=" + "H0*" + prvsig + "[" + (signal_length-2) + "]+H1*" + prvsig + 
                                           "[" + (signal_length-1) + "]+" + "H2*" + prvsig + "[" + 0 + "]+" + "H3*" + prvsig + "[" + 1 + "]; " );
               System.out.println("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
        }
    }

---

Tests
  

A few methods and constants to run some tests.

    public static void test_fwd_cdf44(double[] s, boolean dbg_flag) {
        double[] scopy = new double[s.length];
        System.arraycopy(s, 0, scopy, 0, s.length);
        System.out.print("Input: "); Utils.displaySample(scopy);
        CDF44.orderedDWT(s, dbg_flag);
        System.out.print("FWD CDF(4,4): "); Utils.displaySample(s);
        System.out.println();
    }
   
    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_cdf44(a01a, true);
    }

Here is the output of main. You can plug in different arrays to test it. The value of the forward CDF(4, 4) is the last printed line.

Input: Sample: 4.0 3.0 2.0 1.0
MID           = 2
signal_length = 4
SCL:  s^{1}_{0}[0]=H0*s^{0}_{1}[0]+H1*s^{0}_{1}[1]+H2*s^{0}_{1}[2]+H3*s^{0}_{1}[3]; WVL: s^{1}_{0}[2]=G0*s^{0}_{1}[0]+G1*s^{0}_{1}[1]+G2*s^{0}_{1}[2]+G3*s^{0}_{1}[3]
SCL:  s^{1}_{0}[1]=H0*s^{0}_{1}[2]+H1*s^{0}_{1}[3]+H2*s^{0}_{1}[0]+H3*s^{0}_{1}[1]; WVL: s^{1}_{0}[3]=G0*s^{0}_{1}[2]+G1*s^{0}_{1}[3]+G2*s^{0}_{1}[0]+G3*s^{0}_{1}[1]
FWD CDF(4,4): Sample: 4.760278777324326 2.3107890345411484 -2.220446049250313E-16 1.4142135623730943 

---

Computing & Plotting 3 Scales of the 1D Ordered Haar Wavelet Transform of a Sinusoid Curve: Programmatic Note 2 on Ripples in Mathematics

Vladimir Kulyukin

---

Introduction
  

This is a short programmatic note on Ch. 04, pp. 25 - 27, in "Ripples in Mathematics" by A. Jensen & A. la Cour-Harbo. May the beauty and harmony of this text be conveyed to the readers of this post. These notes are written for those who want to use Java as a programming investigation tool of various wavelet algorithms and Octave to plot the results. I have found this combination of tools to be very useful in my investigations, primarily because both Java and Octave are free. You can find my other notes by searching my blog on "ripples in mathematics." If you are on the same journey, let me know of any bugs in my code or improvements you have found that let us get the results in a better way.

The problem addressed in this post is formulated on p. 25 in "Ripples in Mathematics." The function y = sin(4*pi*t), where t is in [0, 1] is sampled at 512 equidistant points, which gives a discrete signal s9. The index 9, in the formalism adopted in the book, stands for the power of 2 equal to the number of samples in the signal, i.e., 512 = 2^9. This signal is plotted in Figure 4.1 on p. 26. The 1D ordered HWT is applied to this signal over three different scales, i.e., three iterations. The wavelet coefficients are combined into a 1D array in the following order s6, d6, d7, and d8, where d8 is the 1st scale wavelet coefficients, i.e., the wavelet coefficients obtained after the 1st application of the 1D ordered HWT, d7 are the 2nd scale wavelet coefficients (2nd application of the 1D ordered HWT), d6 are the the 3rd scale wavelet coefficients (3rd application of the 1D ordered HWT), and s6 are the averages obtained after the 3rd application of the 1D ordered HWT. This combined array is plotted in Figure 4.2 on p. 26. Figure 4.3 on p. 27 plots s6, d6, d7, d8 individually.

---

Computing the Signal and Wavelet Coefficients in Java


I defined a Java class Ripples_F_p25.java that extends my Function.java class. 

public class Ripples_F_p25 extends Function {
   
    public Ripples_F_p25() {}
    @Override
    public double v(double x) { return Math.sin(4*Math.PI*x/512.0); }
}

Then I defined  RipplesMathCh04.java and used my Partition.java to compute the domain of the sin(4*pi*t/512), apply Ripples_F_p25 to it, and display the range. I defined several static methods (see below): fig_4_1_p26() , fig_4_2_p26() , fig_4_3_d8_p27(),  fig_4_3_d7_p27() , fig_4_3_d6_p27() , and fig_4_3_s6_p27().  These methods compute the ranges for Figures 4.1, 4.2, and 4.3 on pages 26 and 27, respectively, in the book. The methods use several methods of my OneDHaar class to compute the normalized forward and inverse 1D ordered HWTs.
  
public class RipplesInMathCh04 {
  
    public enum SIGNAL { D8, D7, D6, S6 };
  
    static double[] sDomain = Partition.partition(0, 511, 1);
    static double[] sRange  = new double[512];
    static Ripples_F_p25 sRipples_F_p25 = new Ripples_F_p25();
  
    static final int D8_START   = 256;
    static final int D8_END     = 511;
    static final int D7_START   = 128;
    static final int D7_END     = 255;
    static final int D6_START   = 64;
    static final int D6_END     = 127;
    static final int S6_START   = 0; 
    static final int S6_END     = 63;
  
    // prints the range values for the plot in Fig. 4.1, p. 26
    // in "Ripples in Mathematics."
    static void fig_4_1_p26() {
      
        for(int i = 0; i < 512; i++)  {
            sRange[i] = sRipples_F_p25.v(sDomain[i]);
        }
      
        display_signal(sRange);
    }
    // prints the range values for the plot in Fig. 4.2, p.26
    // in "Ripples in Mathematics."
    static void fig_4_2_p26() {
      
        for(int i = 0; i < 512; i++)  {
            sRange[i] = sRipples_F_p25.v(sDomain[i]);
        }
      
        OneDHaar.orderedNormalizedFastHaarWaveletTransformForNumIters(sRange, 3);
      
        for(int i = 0; i < 512; i++)  {
            System.out.println(sRange[i]);
        }
    }
  
    // d8 range values for Fig. 4.3, p. 27 in "Ripples in Mathematics."
    static void fig_4_3_d8_p27() {
      
        for(int i = 0; i < 512; i++)  {
            sRange[i] = sRipples_F_p25.v(sDomain[i]);
        }
      
        OneDHaar.orderedNormalizedFastHaarWaveletTransformForNumIters(sRange, 3);
        display_signal_range(sRange, D8_START, D8_END);
    }
  
    // d7 range values for Fig. 4.3, p. 27 in "Ripples in Mathematics."
    static void fig_4_3_d7_p27() {
      
        for(int i = 0; i < 512; i++)  {
            sRange[i] = sRipples_F_p25.v(sDomain[i]);
        }
      
        OneDHaar.orderedNormalizedFastHaarWaveletTransformForNumIters(sRange, 3);
        display_signal_range(sRange, D7_START, D7_END);
    }
  
    // d6 range values for Fig. 4.3, p. 27 in "Ripples in Mathematics."
    static void fig_4_3_d6_p27() {
      
        for(int i = 0; i < 512; i++)  {
            sRange[i] = sRipples_F_p25.v(sDomain[i]);
        }
      
        OneDHaar.orderedNormalizedFastHaarWaveletTransformForNumIters(sRange, 3);
        display_signal_range(sRange, D6_START, D6_END);
    }
  
    // s6 range values for Fig. 4.3, p. 27 in "Ripples in Mathematics."
    static void fig_4_3_s6_p27() {
      
        for(int i = 0; i < 512; i++)  {
            sRange[i] = sRipples_F_p25.v(sDomain[i]);
        }
      
        OneDHaar.orderedNormalizedFastHaarWaveletTransformForNumIters(sRange, 3);
        display_signal_range(sRange, S6_START, S6_END);
    }}

---

Plotting Java Numbers in Octave


I copied the numbers computed by the above functions from the Java console and pasted them into two Octave scripts ripples_p26.m  and ripples_p27.m that plot the figures on pp 26 and 27 in the book.  The plots produced by these scripts are given below.

---

Integral Approximation with Infinite Limits and Regular Partitions

Vladimir Kulyukin

---

Problem
  

There is a theorem in integral calculus that states that if a function f is integrable on [a, b], then its definite integral on [a, b] is equal to an infinite limit of the sum of f evaluated at regular points obtained with ever decreasing regular partitions. In other words:

Let us investigate how we can use this theorem to programmatically evaluate definite integrals.

---

Function Abstraction

We need to use functions as objects. Toward that result, we define a one-argument Function class that all other subclasses will override. Source code is here.

public class Function {
   
    public Function() {}
    public double v(double x) { return 0; }
    public double[] generateRangeInterval(double[] domain_values) {
        return null;
    }
   
}

The idea of regular partition can now be encapsulated as a subclass of the Function class (RegularPartition.java):

public class RegularPartitionSum extends Function {
   
    private double mA = 0;
    private double mB = 0;
    private Function mF = null;
   
    public RegularPartitionSum(Function f, double a, double b) {
        mA = a;
        mB = b;
        mF = f;
    }
   
    @Override
    public double v(double n) {
        final double step = (mB - mA)/n;
        final int in = (int)n;
        double sum = 0;
        for(int i = 0; i <= in; i++) {
            sum += mF.v(mA + i*step);
        }
        return step*sum;
    }
   
}

Given a function, we should be able to approximate its limit at infinity to any desirable number of steps, which we abstract into a FunctionLimit class (FunctionLimit.java).

public class FunctionLimit {

public static double limitAtInfinity(Function f, double start, double step, int num_steps) {
        double x = start;
        double y = 0;
        while ( num_steps >= 0 ) {
            y = f.v(x);
            x += step;
            num_steps--;
        }
        return y;
}

We proceed by implementing Riemann's sum (RiemannSum.java) with the method that uses FunctionLimit to estimate a definite integral with regular partitions and infinite limits. We compute the limit of the regular partition sums for a specific number of steps.

public class RiemannSum {
    
    public static double infiniteLimitOfRegularPartition(Function f, double a, double b, int num_steps) {
        Function rpF = new RegularPartitionSum(f, a, b);
        return FunctionLimit.limitAtInfinity(rpF, 1, 1, num_steps);
    }
}

---

Sample Functions

Let us define several functions and use the above implementation to evaluate their definite integrals. All of them come from James Stewart's "Calculus: Early Transendentals," 3rd edition. The name of the function class specifies the number of example and the page.

public class F_example_07_p343 extends Function {
   
    public F_example_07_p343() {}
   
    @Override
    public double v(double x) {
        return Math.sqrt(x);
    }
}

The function example_07_p343(), defined below, when placed in the main method, outputs "RS = 4.668915495400732," which is pretty close to the manual evaluation of the definite intergal.
 
public static void example_07_p343() {
        Function f = new F_example_07_p343();
        System.out.println("RS  = " + RiemannSum.infiniteLimitOfRegularPartition(f, 1, 4, 2000));
}

Here is another sample function:

public class F_ex09_p344 extends Function {
   
    public F_ex09_p344() {}
   
    @Override
    public double v(double x) {
        return x*x*x;
    }
}

We can evaluate this function in this method:

public static void ex_09_p344() {
        System.out.println("RS = " + RiemannSum.infiniteLimitOfRegularPartition(new F_ex09_p344(), 0, 5, 1000));
}

The output value is 156.56234375015595. More sample functions are implemented and tested here.

  

---

Midpoint Rule

The accuracy of our evaluations can be compared with the mid-point rule, which is a more standard method of evaluating definite integrals. Add these methods to RiemannSum and use them to evaluate the sample functions.

public static double midPointArea(Function f, double from, double upto, double step) {
        double area = 0;
        if ( from >= upto ) return 0;
        double currP = from;
        while ( currP <= upto - step ) {
            area += f.v(currP+step/2.0);
            currP += step;
        }
        return step*area;
}

public static double midPointRule(Function f, double a, double b, int n) {
        if ( b < a ) throw new IllegalArgumentException("midPointRule: b < a");
        if ( a == b ) return 0;
        final double step = (b - a)/n;
        return RiemannSum.midPointArea(f, a, b, step);
}
  
 

---