Friday, January 15, 2016

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




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