Friday, September 11, 2015

3-Scale 1D Ordered HWT Multiresolution Analysis of a Sinusoid Curve: Programmatic Note 3 on Ripples in Math




Introduction
  
This is my third programmatic note on Ch. 04, pp. 27 - 28, in "Ripples in Mathematics" by A. Jensen & A. la Cour-Harbo. May the beauty and harmony of this text be conveyed to the readers of this post. These notes document my programmatic journey through this text. These notes are written for those who want to use Java as a programming investigation tool of various wavelet algorithms and Octave to plot the results. You can find my other notes by searching my blog on "ripples in mathematics." If you are on the same journey, let me know of any bugs in my code or improvements you have found that let us get the results in a better way.

The problem addressed in this post is formulated on p. 27 in "Ripples in Mathematics." The function y = sin(4*pi*t), where t is in [0, 1] is sampled at 512 equidistant points, which gives a discrete signal s9. Index 9, in the formalism adopted in the book, stands for the power of 2 equal to the number of samples in the signal, i.e., 512 = 2^9

On pp. 28-29, a 3-scale mutiresolution analysis is applied to the sinusoid. In the n-scale multiresolution analysis, we apply n iterations of the 1D ordered HWT and then invert each type of coefficients. The n-scale multiresolution analysis consists of three stages. In the first stage, n iterations of the 1D ordered HWT is applied to the signal. In the second stage, specific coefficients are set to 0's to obtain several arrays. In the third stage, each of the arrays obtained in the second stage is inverted and plotted. 

In the book, the 3-scale multiresolution analysis is applied to s9. The analysis starts on p. 27.  The first stage of the 3-scale multiresolution analysis begins by computing three scales of 1D HWT, which gives us s6, d6, d7, and d8. These coefficients are placed into a 1D array in the following order: s6 first, then d6, d7, and, finally, d8. The second stage of the n-scale multiresolution analysis consists of setting all coefficients in this array with 0's, except for the coefficients that we would like to invert. In other words, first we replace all coefficients in the combined array with 0's except for d8, then we replace all coefficients with 0's except for d7, then all coefficients, except for d6, and, finally, all coefficients, except for s6. In the book, the authors use the notation 06, 07, and 08 to denote the specific regions of the array set to 0. When I was reading the chapter, I decided to modify this notation slightly and used s6 and s06, d6 and d06, d7 and d07, and d8 and d08, where s6, for example, refers to the real s6 coefficients, while s06 refers to all s6 coefficients set to 0. Thus, the array s06, d6, d07, d08 denotes the combined array in which all coefficients are set to 0, except for d6.  This notation helped understand better the explanation on p. 27. Consequently, the second stage of the 3-scale multiresolution analysis gives us four arrays: s06, d06, d07, d8; s06, d06, d7, d08; s06, d6, d07, d08; s6, d06, d07, d08. The third stage of the 3-scale multiresolution analysis consists of inverting each of the four arrays and plotting each one. The result of the 3-scale multiresolution analysis is given in Figure 4.4 on p. 28.

Multiresolution Analysis in Java


I added a method multires_fig_4_4_p28() in RipplesInMathCh04.java that applies the 1D ordered HWT to the signal, sets to 0's all coefficients, except for the values in a specific range, and then prints the result in the Java console. 

public enum SIGNAL { D8, D7, D6, S6 };

static final int D8_START   = 256; 
static final int D8_END       = 511;
static final int D7_START   = 128;
static final int D7_END       = 255;
static final int D6_START   = 64; 
static final int D6_END       = 127;
static final int S6_START   = 0;   
static final int S6_END       = 63;

static double[] sDomain = Partition.partition(0, 511, 1);
static double[] sRange  = new double[512];
static Ripples_F_p25 sRipples_F_p25 = new Ripples_F_p25();

static void multires_fig_4_4_p28(String message, int range_start, int range_end) {
        for(int i = 0; i < 512; i++)  {
            sRange[i] = sRipples_F_p25.v(sDomain[i]);
        }

        OneDHaar.orderedNormalizedFastHaarWaveletTransformForNumIters(sRange, 3);
      
        double[] signal = new double[sRange.length];
        for(int i = 0; i < 512; i++) {
            if ( i >= range_start && i <= range_end ) {
                signal[i] = sRange[i];
            }
            else {
                signal[i] = 0;
            }
        }
      
        System.out.println("=========================");
        System.out.println(message);
        display_signal(signal);
        System.out.println("=========================");
        OneDHaar.orderedNormalizedFastInverseHaarWaveletTransformForNumIters(signal, 3);
        display_signal(signal);
        System.out.println("=========================");
    }


// Fig. 4.4 on p. 28 in "Ripples in Mathematics."
static void fig_4_4_p28(SIGNAL signal) {
        switch ( signal ) {
            case D8: fig_4_4_06_06_07_d8_p28(); break;
            case D7: fig_4_4_06_06_d7_08_p28(); break;
            case D6: fig_4_4_06_d6_07_08_p28(); break;
            case S6: fig_4_4_s6_06_07_08_p28(); break;
        }
    }    


static void fig_4_4_06_06_07_d8_p28() {
        multires_fig_4_4_p28("Fig. 4.4, 06-06-07-d8, p. 28", D8_START, D8_END);
    }
    

static void fig_4_4_06_06_d7_08_p28() {
        multires_fig_4_4_p28("Fig. 4.4, 06-06-d7-08, p. 28", D7_START, D7_END);
    }
    

static void fig_4_4_06_d6_07_08_p28() {
        multires_fig_4_4_p28("Fig. 4.4, 06-d6-07-08, p. 28", D6_START, D6_END);
    }
   

static void fig_4_4_s6_06_07_08_p28() {
        multires_fig_4_4_p28("Fig. 4.4, s6-06-07-08, p. 28", S6_START, S6_END);
    }




Plotting the Results in Octave


I copied the output of the above pasted them into my Octave scripts ripples_s06_d06_d07_d8_p28.m, ripples_s06_d06_d7_d08_p28.m, ripples_s06_d6_d07_d08_p28.m,  and ripples_s6_d06_d07_d08_p28.m. Each of these scripts plots the corresponding 2nd stage array, e.g., s06, d06, d07, d8, and its inverted version given in Figure 4.4 on p. 28. Figure 1 below shows the two graphs plotted by  ripples_s06_d06_d07_d8_p28.m.

Figure 1. Plots produced by  ripples_s06_d06_d07_d8_p28.m