Saturday, October 31, 2015

Relative Mirror Wrapup in DWT: Programmatic Note 6 on Ripples in Mathematics

This is my programmatic note 6 on Ch. 4 in "Ripples in Mathematics" by A. Jensen & A. la Cour-Harbo. As I was implementing CDF(2, 2) defined in Ch 3, I was confronted with what Jensen & Cour-Harbo call the boundary problem. The problem shows up when we need to compute the wavelet values at sub-signal boundaries. This is one of those differences between mathematics and computer science. Mathematicians may have the discrete sample time axis stretching from the negative infinity all the way to the positive infinity. But we, the CS folks, must deal with finite signals of some length n that consist of samples enumerated from 0 up to n-1. How are we, for example, supposed to compute a boundary value if the mathematical equation that we are coding up refers to the samples that do not exist.  Say, computing the value at sample 0 may require the values at samples -1 and -2. This is the left boundary problem. Similarly, suppose the length of our signal is n and we need to compute the value at sample n-1 for which the equation requires us to access the samples at positions n and n+1, which do not exist. 

One common solution is to return 0's for the samples outside our sub-signal's range. Professor Nievergelt in his excellent book "Wavelets Made Easy" shows that this may not work as well as we would naively expect. Jensen & Cour-Harbo do not seem to explain which method they used in their CDF(2, 2) experiments in Ch. 4. One other method suggested in "Wavelets Made Easy" is mirror reflection. This is what I decided to implement. I may change my implementation in the future as I run experiments with large data sets. But, for now, I will keep experimenting with what I call relative mirror wrapup which I briefly explain below.

Relative Mirror Wrapup
 Suppose we have a signal i.e., a 1D array of values: [a(0), a(1), a(2), a(3), a(4), a(5), a(6), a(7)]. When we compute some DWT, especially its ordered version, we must divide the signal into sub-signals and apply the DWT to each sub-signal according to some recurrence. So, we have the start and end of our sub-signal, e.g., s and e, and then we have the discrete sample axis stretching from the negative infinite to the positive infinity. We can align our sub-signal's start (s) with the 0 of the sample axis. Hence, the term relative, because the 0 on the discrete sample axis is positioned relative to the start of the sub-signal. For example,

sample axis:        -3  -2  -1   0       1      2       3      4      5       6      7     8 ....  
signal:                              [a(0), a(1), a(2), a(3), a(4), a(5), a(6), a(7)]
sub-signal's start & end:    s       e

The above notation states that our sub-signal is [a(0), a(1)]. Its starting index, s, is 0 and its end index, e, is 1. The notation a(k), 0 <= k <= 7, denotes the individual samples in our signal where k is an array index. The sample axis above the signal is placed so that the start of the sub-signal is aligned with its 0. 

Let us define a method that allows us to map the values from the sample axis into the sub-signal indices using the mirror principle. We start with the left boundary. Suppose we need to access a(-1). If we put the mirror at the left boundary of the sub-signal, the reflection of a(-1) is a(0). In other words, the sample axis value of -1 is mapped to the sub-signal's index 0. What about computing a(-3)? The mirror at the left boundary reflects -2 as 1. Similarly, -3 is reflected to 0, and so on. 

If we place the mirror at the right boundary, we reflect a(2) to a(1). Note that we cannot simply return a(2), although it really exists in the larger signal, because our sub-signal ends at 1. Consequently, all indices on the sample axis less than 0 and greater than 1 must be reflected back into the sub-signal's indices. Let us consider another example.

sample axis:      -3   -2   -1    0       1       2      3      4       5       6      7     8 ....  
signal:                               [a(0), a(1), a(2), a(3), a(4), a(5), a(6), a(7)]
sub-signal's start & end:              s                        e
The beginning of the signal is aligned with the sample axis. The sub-signal starts at 1 and ends at 4. If we place the mirror at the left boundary of the sub-signal, i.e., position 1, then a(0) is reflected to a(1), a(-1) to a(2), a(-2) to a(3), a(-3) to a(4) and then the reflection pattern repeats. If the mirror is placed at the right boundary of the sub-signal, a(5) is reflected to a(4), a(6) to a(3), a(7) to a(2), a(8) to a(1) and then the reflection pattern repeats.

Java Implementation
I have implemented the above logic in the Java class given below.

public class DWTBoundaryHandler {
    // relative mirror wrapup
    public static int rmw(int s, int e, int i) {
        if ( s > e ) throw new IllegalArgumentException("Start = " + s + " > End = " + e);
        final int n = (e - s) + 1;
        if ( i >= 0 && i < n ) {
            return s+i;
        else {
            int r = Math.abs(i) % n;
            if ( i < 0 ) {
                return ( r != 0 )? s + r -1: e;
            else {
                return ( r != 0 )? e - r: e;
    static void run_rmw_test(int s, int e, int is, int ie) { for(int i = is; i <= ie; i++) { display_rmw(s, e, i); }}
    // -3  -2  -1    0      1      2       3      4      5       6      7     8 .... 
    //              [a(0), a(1), a(2), a(3), a(4), a(5), a(6), a(7)]
    //                  s      e

    public static void rmw_test00() { run_rmw_test(0, 1, -4, -1); run_rmw_test(0, 1,  2,  5); }
    // -3  -2  -1    0       1       2      3       4      5      6       7    8 .... 
    //               [a(0), a(1), a(2), a(3), a(4), a(5), a(6), a(7)]
    //                          s=e

    public static void rmw_test01(int s) {
        display_rmw(s, s, -1);
        display_rmw(s, s, -2);
        display_rmw(s, s, -3);
        display_rmw(s, s, 0);
        display_rmw(s, s, 1);
        display_rmw(s, s, 2);
        display_rmw(s, s, 3);
        display_rmw(s, s, 4);
    //    -3   -2   -1     0      1      2      3       4     5     6     7    8 .... 
    //                    [a(0), a(1), a(2), a(3), a(4), a(5), a(6), a(7)]
    //                                 s                        e

    public static void rmw_test02() { run_rmw_test(1, 4, -5, -1); run_rmw_test(1, 4, 4,  8); }
    public static void display_rmw(int s, int e, int i) {
        System.out.println("rmw(s=" + s + ",e=" + e + ",i=" + i + ")=" + DWTBoundaryHandler.rmw(s, e, i));
    public static void main(String[] args) { rmw_test02(); }