VKEDCO: Vladimir Kulyukin's Education Coop: cdf

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

Vladimir Kulyukin

Introduction

This is a conceptual note on Ch. 3, pp. 22-23, in "Ripples in Mathematics" by A. Jensen & A. la Cour-Harbo. I had implemented the 1D and 2D Haar Transforms (my github repo is here) before reading this great and beautiful book. So, on my first reading, I quickly read through the notation on p. 22. However, when I started implementing CDF(2, 2), the issue of notation became important. Below are the modifications that I introduced to the DWT scale notation adopted in Ch 3 to understand CDF(2, 2). The notation for the forward DWT used in Ch. 3 is given in Fig. 1.

Figure 1. Notation for Forward DWT

The forward direct wavelet transform (DWT) is applied to a signal of length 2^J, for J > 0. Each application of the transform is called a scale. The maximum scale value is J. I take scale 0 to be the original signal. A signal is a sequence of samples, e.g., real numbers. At each scale greater than 0, the signal can be represented as the signal proper, i.e., the scaled signal, and the difference sequence (this term is mine). The difference sequence consists of numbers that allow us to recover the signal at the previous scale and to reason about the signal properties, e.g., how the signal changes from the 1st half to the second half. Fig. 2 shows the notation I used to represent a signal at a specific scale to make sense of the CDF(2, 2) definition on p. 23. This notation differs slightly from the one used in the book and may be helpful to the readers of this document when they implement CDF(2, 2).

Figure 2. Notation for Signal at Scale L

Fig. 3 gives my notation for the i-th sample of a scaled signal. In other words, the symbol defined in Fig. 2 denotes a signal whose individual elements are referenced by symbols defined in Fig. 3. Note that in Figures 3 and 5, the range for i should be defined as 0 <= i < 2^(J-L), not 0 <= i <= 2^(J-L). I noticed this typo after I published this post.

Figure 3. Notation for the I-th Sample of the Signal at Scale L

The second component of the DWT at a given scale is the difference sequence. I do not think difference sequence is a standard term. The notation for difference sequences is shown in Fig. 4.