VKEDCO: Vladimir Kulyukin's Education Coop: January 2016

Wednesday, January 20, 2016

Computing Inverse Ordered CDF(4, 4): Programmatic Note 9 on Ripples in Mathematics

Problem

This is my programmatic note 9 on "Ripples in Mathematics" by A. Jensen & A. la Cour-Harbo. In this note, I detail my Java implementation of the inverse ordered CDF(4,4) or D4 as it is known in the wavelet literature. The ordered forward CDF(4,4) is covered in programmatic note 7. My Java source is here.

Inverse Ordered CDF(4, 4) in Java

I will skip the definition of the CDF44 class for the sake of space and give only the definition of the static method that computes inverse ordered D4. The method is a static method of CDF44. Utils.isPowerOf2(N) is from my Utils class. It has a bunch of util methods, primarily numerical.

     public static void orderedInverseDWT(double[] signal_transform, boolean dbg_flag) {
        final int N = signal_transform.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 numInvScalesToDo = Utils.powVal(N)-1;
        int currInvScale     = 0;
        int transform_length = 4;

        if ( dbg_flag ) {
            System.out.print("<=INPUT: "); Utils.displaySample(signal_transform);
        }

        while ( transform_length <= N ) {
            int mid = transform_length >> 1;
            if ( dbg_flag ) System.out.println("MID              = " + mid);
            if ( dbg_flag ) System.out.println("transform_length = " + transform_length);

            double[] inv_sig = new double[transform_length]; // restored values

            String cur_sig = null;
            String prv_sig = null;

            if ( dbg_flag ) {
                cur_sig = "s^{" + (numInvScalesToDo-1) + "}_{" + (currInvScale+1) + "}";
                prv_sig = "s^{" + numInvScalesToDo     + "}_{" + currInvScale     + "}";
            }

            inv_sig[0] = IH0*signal_transform[mid-1] + IH1*signal_transform[transform_length-1] + IH2*signal_transform[0] + IH3*signal_transform[mid];
            inv_sig[1] = IG0*signal_transform[mid-1] + IG1*signal_transform[transform_length-1] + IG2*signal_transform[0] + IG3*signal_transform[mid];

            if ( dbg_flag ) {
                System.out.println("INV SCL: " + cur_sig + "[" + 0 + "] = " + "IH0*" + prv_sig + "[" + (mid-1) + "] + " +
                        "IH1*" + prv_sig + "[" + (transform_length-1) + "] + " + "IH2*" + prv_sig + "[0] + " + "IH3*" + prv_sig + "[" + mid + "]");
                System.out.println("INV WVL: " + cur_sig + "[" + 1 + "] = " + "IG0*" + prv_sig + "[" + (mid-1) + "] + " +
                        "IG1*" + prv_sig + "[" + (transform_length-1) + "] + " + "IG2*" + prv_sig + "[0] + " + "IG3*" + prv_sig + "[" + mid + "]");
            }


            int i = 0, j = 2;

            while ( i < mid-1 ) {
                if ( dbg_flag ) {
                    cur_sig = "s^{" + (numInvScalesToDo-1) + "}_{" + (currInvScale+1) + "}";
                    prv_sig = "s^{" + numInvScalesToDo     + "}_{" + currInvScale     + "}";
                }

                if ( dbg_flag ) {
                    System.out.println("INV SCL: " + cur_sig + "[" + j + "] = " + "IH0*" + prv_sig + "[" + i + "] + " +
                        "IH1*" + prv_sig + "[" + (mid+i) + "] + " + "IH2*" + prv_sig + "[" + (i+1) + "] + " +
                        "IH3*" + prv_sig + "[" + (mid+i+1) + "]");
                }

                //           scalers                     wavelets
                inv_sig[j] = IH0*signal_transform[i]   + IH1*signal_transform[mid+i] +
                         IH2*signal_transform[i+1] + IH3*signal_transform[mid+i+1];

                if ( dbg_flag ) {
                    System.out.println("INV WVL: " + cur_sig + "[" + (j+1) + "] = " + "IG0*" + prv_sig + "[" + i + "] + " +
                        "IG1*" + prv_sig + "[" + (mid+i) + "] + " + "IG2*" + prv_sig + "[" + (i+1) + "] + " +
                        "IG3*" + prv_sig + "[" + (mid+i+1) + "]");
                }

                //             scalers                     wavelets
                inv_sig[j+1] = IG0*signal_transform[i]   + IG1*signal_transform[mid+i] +
                               IG2*signal_transform[i+1] + IG3*signal_transform[mid+i+1];

                i += 1; j += 2;
            }

            currInvScale     += 1;
            numInvScalesToDo -= 1;

            System.arraycopy(inv_sig, 0, signal_transform, 0, inv_sig.length);
            transform_length <<= 1; // multiply by 2
        }
    }

Tests

A few methods and constants to run some tests.The method test_fwd_inv_cdf44 runs the ordered forward D4 on the signal s by calling CDF44.orderedDWT(s, dbg_flag) and then inverses the transform by calling CDF44.orderedInverseDWT(s, dbg_flag).

    public static void test_fwd_inv_cdf44(double[] s, boolean dbg_flag) {
        System.out.print("Input: "); Utils.displaySample(s);

        CDF44.orderedDWT(s, dbg_flag);

        System.out.print("FWD CDF(4,4): "); Utils.displaySample(s);
        System.out.println();

        CDF44.orderedInverseDWT(s, dbg_flag);

        System.out.print("INV 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_inv_cdf44(a02a, false);
        test_fwd_inv_cdf44(a02b, false);
    }

Here is the output of main. You can plug in different arrays to test the transform. The debug flag is set to false.

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

Friday, January 15, 2016

Computing Forward Ordered CDF(4, 4) for a Given Number of Iterations (Scales): Programmatic Note 8 on Ripples in Mathematics

Vladimir Kulyukin

Problem

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.

Forward CDF(4, 4) for a Given Number of Iterations/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*" +

                                                                          prvsig + "[" + (signal_length-1) + "]+" +
                                                                          "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;
        }
    }

Tests

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

Thursday, January 7, 2016

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