Thursday, January 7, 2016

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




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