Monday, September 14, 2015

Noise Detection in 3-Scale Multiresolution HWT Analysis of a Synthetic Sinusoid: Programmatic Note 4 on Ripples in Mathematics

This is my programmatic note 4 on Ch. 04, p. 29, 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.

The problem addressed in this post is formulated on p. 28 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. Then some noise is introduced into the signal by setting the value of the 200-th sample to 2. The question is, at which scale of the 1D ordered HWT can the noise be detected?

The problem is solved by applying three iterations of the 1D ordered HWT to the modified signal. Then plotting each of the components of the resultant signal, i.e.., d8, d7, d6, and s6 and observing which plot gives us a spike at sample 200. This is what is done in Figure 4.6 on p. 29 in the book.

Computing 3 Iterations of 1D ordered HWT in Java

The methods fig_4_6_s06_d06_d07_d8_p29(), fig_4_6_s06_d06_d7_d08_p29(),  fig_4_6_s06_d6_d07_d08_p29(), fig_4_6_s06_d6_d07_d08_p29(), fig_4_6_s6_d06_d07_d08_p29()  in 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.

static void fig_4_6_s06_d06_d07_d8_p29() { multires_fig_4_6_p29("06-06-07-d8, Fig. 4.6, p. 29", D8_START, D8_END); }
static void fig_4_6_s06_d06_d7_d08_p29() { multires_fig_4_6_p29("06-06-d7-08, Fig. 4.6, p. 29", D7_START, D7_END); }
static void fig_4_6_s06_d6_d07_d08_p29() { multires_fig_4_6_p29("06-d6-07-08, Fig. 4.6, p. 29", D6_START, D6_END); }
static void fig_4_6_s6_d06_d07_d08_p29() { multires_fig_4_6_p29("s6-06-07-08, Fig. 4.6, p. 29", S6_START, S6_END); }

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

        sRange[200] = 2; // spike at 200
        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;
        OneDHaar.orderedNormalizedFastInverseHaarWaveletTransformForNumIters(signal, 3);


Plotting D6, D7, D6, & S6 in Octave

As in my previous programmatic notes on "Ripples," I took the numbers computed by Java and copied and pasted them into several Octave scripts : ripples_s06_d06_d07_d8.m, ripples_s06_d06_d7_d08.m , ripples_s06_d6_d07_d08.m , ripples_s6_d06_d07_d08.m , and ripples_fig_4_6_p29.m . You can refer to this post for a brief explanation of my notation such as s06, d06, d07, etc., which is slightly different from the one adopted in the book. Figures 1 and 2 below give the outputs of two scripts.

Fig. 1. Plots produced by ripples_s06_d6_d07_d08.m

Fig. 2. Plots produced by ripples_fig_4_5_p29.m
As can be seen from Fig. 2, the noise is easily detectable in D8, D7, or D6. It is barely noticeable in S6. So it does make a difference at which scale noise detection is done.