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