## Friday, August 28, 2015

### Reconstruction of Thresholded 1D Signals Processed with 1D Ordered Haar Wavelet Transform: Programmatic Note 1 on Ripples in Math

Introduction

This is a short programmatic note on Ch. 02, pp. 7 - 9, in "Ripples in Mathematics" by A. Jensen & A. la Cour-Harbo. I am writing it with a deep gratitude to the authors for such a wonderful text. I have read the text from cover to cover and enjoyed every page. Although I am still fuzzy on some technical details, I hope to gain some understanding as I re-read and programmatically re-construct some concepts covered in this great book. May the beauty and harmony of this text be conveyed to the readers of this post. 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 described in this post is formulated on p. 7. The problem is to compute 3 forward iterations of the 1D ordered Haar transform of an 8-sample signal [56, 40, 8, 24, 48, 48, 40, 16].

Computing Table 2.1, p. 7, Table 2.2, p. 8, and Table 2.3, p. 9

Table 2.1 on p. 7 gives the above 8-element signal and 3 consecutive applications of the 1D Haar wavelet to it. Here it is.

------------------------------------------------------------------------
56.00    40.00    8.00    24.00    48.00    48.00    40.00    16.00
---------------------------------------------------------------------------
48.00    16.00    48.00    28.00    8.00    -8.00    0.00    12.00
---------------------------------------------------------------------------
32.00    38.00    16.00    10.00    8.00    -8.00    0.00    12.00
---------------------------------------------------------------------------
35.00    -3.00    16.00    10.00    8.00    -8.00    0.00    12.00
---------------------------------------------------------------------------

I have written a Java class, RipplesInMathCh02.java, to construct this table. The class uses my implementation of the 1D Haar filter available at here.

static void computeTableForSignalOnPage7(double[] signal, double thresh, String tableTitle) {
System.out.println(tableTitle);

final int n = signal.length;
double[] signalFor1InverseIter  = new double[n];
double[] signalFor2InverseIters = new double[n];
double[] signalFor3InverseIters = new double[n];
double[] signalFor3Iters             = new double[n];

System.arraycopy(signal, 0, signalFor3Iters, 0, n);
OneDHaar.orderedFastHaarWaveletTransformForNumIters(signalFor3Iters, 3);

System.arraycopy(signalFor3Iters, 0, signalFor3InverseIters,  0, n);
OneDHaar.orderedFastInverseHaarWaveletTransformForNumIters(signalFor3InverseIters, 3, thresh);

System.arraycopy(signalFor3Iters, 0, signalFor2InverseIters,  0, n);
OneDHaar.orderedFastInverseHaarWaveletTransformForNumIters(signalFor2InverseIters, 2, thresh);

System.arraycopy(signalFor3Iters, 0, signalFor1InverseIter,  0, n);
OneDHaar.orderedFastInverseHaarWaveletTransformForNumIters(signalFor1InverseIter, 1, thresh);

OneDHaar.thresholdSignal(signalFor3Iters, thresh);

printTable(signalFor3InverseIters, signalFor2InverseIters, signalFor1InverseIter, signalFor3Iters, 62);

}

The array signalFor3Iters is the original 8-element signal to which we apply three forward iterations of the 1D ordered Fast Haar Wavelet Transform (FHWT). I use System.arraycopy, because my implementation of the 1D ordered FHWT is destructive, i.e., it destructively modifies the original signal.

System.arraycopy(signal, 0, signalFor3Iters, 0, n);
OneDHaar.orderedFastHaarWaveletTransformForNumIters(signalFor3Iters, 3);

After the FHWT is applied, I apply the inverse FHWT for a specific number of iterations and a specific threshold. For example, OneDHaar.orderedFastInverseHaarWaveletTransformForNumIters(signalFor3InverseIters, 3, thresh) applies 3 iterations of the inverse FHWT to signalFor3InverseIters where all samples whose absolute value is less then or equal to the threshold are set to 0. In other words,

public static void thresholdSignal(double[] signal, double thresh) {
final int n = signal.length;
double thresholdedSignal[] = new double[n];
for(int t = 0; t < n; t++) {
if ( Math.abs(signal[t]) > thresh )
thresholdedSignal[t] = signal[t];
else
thresholdedSignal[t] = 0;
}

System.arraycopy(thresholdedSignal, 0, signal, 0, n);
thresholdedSignal = null;

}

The method  prinTable, which is the last call in the computeTableForSignalOnPage7() method, simply formats and displays the computational results in the format that matches the one in the text. I have also implemented the version of computeTableForSignalOnPage7() that does not threshold the signal.

static void computeTableForSignalOnPage7(double[] signal, String tableTitle) {
System.out.println(tableTitle);

final int n = signal.length;
double[] signalFor1InverseIter  = new double[n];
double[] signalFor2InverseIters = new double[n];
double[] signalFor3InverseIters = new double[n];
double[] signalFor3Iters        = new double[n];

System.arraycopy(signal, 0, signalFor3Iters, 0, n);
OneDHaar.orderedFastHaarWaveletTransformForNumIters(signalFor3Iters, 3);

System.arraycopy(signalFor3Iters, 0, signalFor3InverseIters,  0, n);
OneDHaar.orderedFastInverseHaarWaveletTransformForNumIters(signalFor3InverseIters, 3);

System.arraycopy(signalFor3Iters, 0, signalFor2InverseIters,  0, n);
OneDHaar.orderedFastInverseHaarWaveletTransformForNumIters(signalFor2InverseIters, 2);

System.arraycopy(signalFor3Iters, 0, signalFor1InverseIter,  0, n);
OneDHaar.orderedFastInverseHaarWaveletTransformForNumIters(signalFor1InverseIter, 1);

printTable(signalFor3InverseIters, signalFor2InverseIters, signalFor1InverseIter, signalFor3Iters, 62);

}

With these methods in place, I implemented the methods that compute Tables 2.1., 2.2. and 2.3 on pp. 7 - 9.

static void computeTable_2_1_p7(double[] signal) {
computeTableForSignalOnPage7(signal, "Table 2.1, p. 7, Ch. 2, \"Ripples in Mathematics\"" );

}
static void computeTable_2_2_p8(double[] signal) {
computeTableForSignalOnPage7(signal, 4, "Table 2.2, p. 7, Ch. 8, \"Ripples in Mathematics\"" );

}
static void computeTable_2_3_p9(double[] signal) {
computeTableForSignalOnPage7(signal, 9, "Table 2.3, p. 9, Ch. 2, \"Ripples in Mathematics\"" );

}

I wrapped these 3 methods for clarity into the methods that explicitly specify the page numbers and then placed them in the main():

static void page7() { computeTable_2_1_p7(signal_p7); }
static void page8() { computeTable_2_2_p8(signal_p7); }
static void page9() { computeTable_2_3_p9(signal_p7); }

public static void main(String[] args) {
page7();
page8();
page9();

}

The output of the main() above is as follows:

Table 2.1, p. 7, Ch. 2, "Ripples in Mathematics"
--------------------------------------------------------------
56.00    40.00    8.00    24.00    48.00    48.00    40.00    16.00
--------------------------------------------------------------
48.00    16.00    48.00    28.00    8.00    -8.00    0.00    12.00
--------------------------------------------------------------
32.00    38.00    16.00    10.00    8.00    -8.00    0.00    12.00
--------------------------------------------------------------
35.00    -3.00    16.00    10.00    8.00    -8.00    0.00    12.00
--------------------------------------------------------------

Table 2.2, p. 7, Ch. 2, "Ripples in Mathematics"
--------------------------------------------------------------
59.00    43.00    11.00    27.00    45.00    45.00    37.00    13.00
--------------------------------------------------------------
51.00    19.00    45.00    25.00    8.00    -8.00    0.00    12.00
--------------------------------------------------------------
35.00    35.00    16.00    10.00    8.00    -8.00    0.00    12.00
--------------------------------------------------------------
35.00    0.00    16.00    10.00    8.00    -8.00    0.00    12.00
--------------------------------------------------------------

Table 2.2, p. 7, Ch. 2, "Ripples in Mathematics"
--------------------------------------------------------------
51.00    51.00    19.00    19.00    45.00    45.00    37.00    13.00
--------------------------------------------------------------
51.00    19.00    45.00    25.00    0.00    0.00    0.00    12.00
--------------------------------------------------------------
35.00    35.00    16.00    10.00    0.00    0.00    0.00    12.00
--------------------------------------------------------------
35.00    0.00    16.00    10.00    0.00    0.00    0.00    12.00
--------------------------------------------------------------

Figures 2.1 & 2.2 on p. 9

I placed the above numbers into an Octave script, ripples_p9.m, and graphed them. When I run this script, I get the following plots which match the ones given in the text.

## Friday, August 14, 2015

### 3D Plots of Toroidal Spirals in Octave

Problem

A toroidal spiral is a space curve whose x, y, and z coordinates are computed with the following parametric equations:
x = (C + sin(Kt))cos(t); y = (C + sin(Kt))sin(t); z = cos(20t). They are called toroidal because they lie on a torus, as shown in Figure 1. In this post, we will briefly tackle the problem of 3D plotting toroidal spirals in Octave and investigate the impact of the C and K parameters.

Figure 1. 3D plot of a Toroidal Spiral

Solution Outline

The Octave code below shows how the plot in Figure 1 can be generated in Octave. It also shows a generic approach and allows us to play with the C and K parameters.

t = 0:0.01:2*pi;
K = 20;
C = 2;

x2 = arrayfun(@(x) (C + sin(K*x))*cos(x), t);
y2 = arrayfun(@(x) (C + sin(K*x))*sin(x), t);
z2 = arrayfun(@(x) cos(K*x), t);

figure;
plot3(x2, y2, z2);
xlabel('x=(2+sin(20t))cos(t)');
ylabel('y=(2+sin(20t))sin(t)');
zlabel('z=cos(20t)');
view([30, 40]);

Figure 2 shows another 3D plot  with C=4 and K=20.

Figure 2. 3D Plot of Another Toroidal Curve

Octave Source

t = 0:0.01:2*pi;
K = 20;

C = 0;
x0 = arrayfun(@(x) (C + sin(K*x))*cos(x), t);
y0 = arrayfun(@(x) (C + sin(K*x))*sin(x), t);
z0 = arrayfun(@(x) cos(K*x), t);

figure;
plot3(x0, y0, z0);
xlabel('x=(0+sin(20t))cos(t)');
ylabel('y=(0+sin(20t))sin(t)');
zlabel('z=cos(20t)');
view([30, 40]);

C = 1;
x1 = arrayfun(@(x) (C + sin(K*x))*cos(x), t);
y1 = arrayfun(@(x) (C + sin(K*x))*sin(x), t);
z1 = arrayfun(@(x) cos(K*x), t);

figure;
plot3(x1, y1, z1);
xlabel('x=(1+sin(20t))cos(t)');
ylabel('y=(1+sin(20t))sin(t)');
zlabel('z=cos(20t)');
view([30, 40]);

C = 2;
x2 = arrayfun(@(x) (C + sin(K*x))*cos(x), t);
y2 = arrayfun(@(x) (C + sin(K*x))*sin(x), t);
z2 = arrayfun(@(x) cos(K*x), t);

figure;
plot3(x2, y2, z2);
xlabel('x=(2+sin(20t))cos(t)');
ylabel('y=(2+sin(20t))sin(t)');
zlabel('z=cos(20t)');
view([30, 40]);

C = 3;
x3 = arrayfun(@(x) (C + sin(K*x))*cos(x), t);
y3 = arrayfun(@(x) (C + sin(K*x))*sin(x), t);
z3 = arrayfun(@(x) cos(K*x), t);

figure;
plot3(x3, y3, z3);
xlabel('x=(3+sin(20t))cos(t)');
ylabel('y=(3+sin(20t))sin(t)');
zlabel('z=cos(20t)');
view([30, 40]);

C = 4;
x4 = arrayfun(@(x) (C + sin(K*x))*cos(x), t);
y4 = arrayfun(@(x) (C + sin(K*x))*sin(x), t);
z4 = arrayfun(@(x) cos(K*x), t);

figure;
plot3(x4, y4, z4);
xlabel('x=(4+sin(20t))cos(t)');
ylabel('y=(4+sin(20t))sin(t)');
zlabel('z=cos(20t)');
view([30, 40]);

C = 5;
x5 = arrayfun(@(x) (CONST + sin(K*x))*cos(x), t);
y5 = arrayfun(@(x) (CONST + sin(K*x))*sin(x), t);
z5 = arrayfun(@(x) cos(K*x), t);

figure;
plot3(x5, y5, z5);
xlabel('x=(5+sin(20t))cos(t)');
ylabel('y=(5+sin(20t))sin(t)');
zlabel('z=cos(20t)');
view([30, 40]);