Demonstrate how the Haar wavelet transform can be calculated using matrices.
The even rows of a forward transform matrix calculate the Haar wavelet scaling function. This is a smoothed value, which in the case of the Haar transform is simply the average of two values. The scaling function is sometimes referred to as a low pass filter.
The odd rows in the forward transform matrix calculate the Haar wavelet function. This represents the change between two elements. In the case of this version of the Haar algorithm, this is the average difference between two adjacent values. The wavelet function is sometimes referred to as a high pass filter.
*/ class haarDemo { // // Haar wavelet forward transform matrices // static final double fwdMat16[][] = { {0.5, 0.5, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}, {0.5,-0.5, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}, {0, 0, 0.5, 0.5, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}, {0, 0, 0.5,-0.5, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}, {0, 0, 0, 0, 0.5, 0.5, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}, {0, 0, 0, 0, 0.5,-0.5, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}, {0, 0, 0, 0, 0, 0, 0.5, 0.5, 0, 0, 0, 0, 0, 0, 0, 0}, {0, 0, 0, 0, 0, 0, 0.5,-0.5, 0, 0, 0, 0, 0, 0, 0, 0}, {0, 0, 0, 0, 0, 0, 0, 0, 0.5, 0.5, 0, 0, 0, 0, 0, 0}, {0, 0, 0, 0, 0, 0, 0, 0, 0.5,-0.5, 0, 0, 0, 0, 0, 0}, {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0.5, 0.5, 0, 0, 0, 0}, {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0.5,-0.5, 0, 0, 0, 0}, {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0.5, 0.5, 0, 0}, {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0.5,-0.5, 0, 0}, {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0.5, 0.5}, {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0.5,-0.5}}; static final double fwdMat8[][] = {{0.5, 0.5, 0, 0, 0, 0, 0, 0}, {0.5,-0.5, 0, 0, 0, 0, 0, 0}, {0, 0, 0.5, 0.5, 0, 0, 0, 0}, {0, 0, 0.5,-0.5, 0, 0, 0, 0}, {0, 0, 0, 0, 0.5, 0.5, 0, 0}, {0, 0, 0, 0, 0.5,-0.5, 0, 0}, {0, 0, 0, 0, 0, 0, 0.5, 0.5}, {0, 0, 0, 0, 0, 0, 0.5,-0.5}}; static final double fwdMat4[][] = {{0.5, 0.5, 0, 0}, {0.5,-0.5, 0, 0}, {0, 0, 0.5, 0.5}, {0, 0, 0.5,-0.5}}; static final double fwdMat2[][] = {{0.5, 0.5}, {0.5,-0.5}}; // // Haar inverse transform matrices // static final double invMat16[][] = { {1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}, {1,-1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}, {0, 0, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}, {0, 0, 1,-1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}, {0, 0, 0, 0, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}, {0, 0, 0, 0, 1,-1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}, {0, 0, 0, 0, 0, 0, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0}, {0, 0, 0, 0, 0, 0, 1,-1, 0, 0, 0, 0, 0, 0, 0, 0}, {0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 0, 0, 0, 0, 0, 0}, {0, 0, 0, 0, 0, 0, 0, 0, 1,-1, 0, 0, 0, 0, 0, 0}, {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 0, 0, 0, 0}, {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1,-1, 0, 0, 0, 0}, {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 0, 0}, {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1,-1, 0, 0}, {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1}, {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1,-1}}; static final double invMat8[][] = {{1, 1, 0, 0, 0, 0, 0, 0}, {1,-1, 0, 0, 0, 0, 0, 0}, {0, 0, 1, 1, 0, 0, 0, 0}, {0, 0, 1,-1, 0, 0, 0, 0}, {0, 0, 0, 0, 1, 1, 0, 0}, {0, 0, 0, 0, 1,-1, 0, 0}, {0, 0, 0, 0, 0, 0, 1, 1}, {0, 0, 0, 0, 0, 0, 1,-1}}; static final double invMat4[][] = {{1, 1, 0, 0}, {1,-1, 0, 0}, {0, 0, 1, 1}, {0, 0, 1,-1}}; static final double invMat2[][] = {{1, 1}, {1,-1}}; public static void main( String argv[] ) { double v[] = {32.0, 10.0, 20.0, 38.0, 37.0, 28.0, 38.0, 34.0, 8.0, 24.0, 18.0, 9.0, 23.0, 24.0, 28.0, 34.0}; matWavelet w = new matWavelet(); linearPr p = new linearPr(); System.out.println("before:"); p.prVec( v ); System.out.println(); w.transStep(fwdMat16, v ); w.transStep(fwdMat8, v ); w.transStep(fwdMat4, v ); w.transStep(fwdMat2, v ); p.pr_ordered( v ); System.out.println(); w.invTransStep( invMat2, v ); w.invTransStep( invMat4, v ); w.invTransStep( invMat8, v ); w.invTransStep( invMat16, v ); System.out.println("after:"); p.prVec( v ); } // main } // haarDemo matrix/linear/ 40755 1040 1001 0 7422407574 11745 5 ustar everyone matrix/linear/linearOps.java 100644 1040 1001 17450 7422401351 14653 0 ustar everyone package linear; import java.lang.Math.*; /**Linear algebra (vector, vector/matrix and matrix) operations
The linear algebra operations are implemented for vectors and matrices of double elements.
Many of the functions are overloaded. In one case a new array or matrix is allocated for the function result. In the other case the result is stored in an argument.
This is not an "industrial" quality linear algebra package, since the functions do not do enough checking for array arguments with null values or length mismatch.
author: Ian Kaplan
use: You may use this software for any purpose as
long as the author, Ian Kaplan, cannot be held liable for
the result. Please include the author attribution.
Calculate the inner, or dot, product of two vectors. The result of a dot product operation is a double value. To produce a correct result, the two vectors should be the same length. To avoid indexing errors, this function will calculate a result using the minimum vector length.
*/ public double dotprod( double vecA[], double vecB[] ) { int len = Math.min( vecA.length, vecB.length ); double innerProd = 0.0; for (int i = 0; i < len; i++) { innerProd = innerProd + vecA[i] * vecB[i]; } return innerProd; } // dotprod /**Vector magnitude: ||v||
Vector magnitude is the square of a vector, which is the same as the dot product of the vector with itself.
*/ public double vecmag( double v[] ) { return dotprod( v, v ); } // vecmag /**Add two vectors, a and b. The result is returned in a new vector. The two vectors should be the same length. If this is not the case, a result that corresponds to the smallest vector will be returned.
*/ public double[] vecadd( double a[], double b[] ) { int len = Math.min(a.length, b.length ); double rslt[] = new double[len]; for (int i = 0; i < len; i++) { rslt[i] = a[i] + b[i]; } return rslt; } // vecadd /**Add two vectors, a and b. The result is returned in argument c.
*/ public void vecadd( double c[], double a[], double b[] ) { int len = Math.min(a.length, b.length ); for (int i = 0; i < len; i++) { c[i] = a[i] + b[i]; } } // vecadd /**Subtract vector b from vector a. The result is returned in a new vector.
*/ public double[] vecsub( double a[], double b[] ) { int len = Math.min(a.length, b.length ); double rslt[] = new double[len]; for (int i = 0; i < len; i++) { rslt[i] = a[i] - b[i]; } return rslt; } // vecsub /**Subtract vector b from vector a. The result is returned in the argument c.
*/ public void vecsub(double c[], double a[], double b[] ) { int len = Math.min(a.length, b.length ); for (int i = 0; i < len; i++) { c[i] = a[i] - b[i]; } } // vecsub /** Multiply a vector, v, by a scalar value, scale. The result is returned in a new vector. */ public double[] multScaleVec( double v[], double scale ) { int len = v.length; double rslt[] = new double[len]; for (int i = 0; i < len; i++) { rslt[i] = v[i] * scale; } return rslt; } // multScaleVec /**Matrix multiply
This is the linear algebra text book function to mltiply two matrices, where matA is MxN and matB is NxM, yielding an MxM result. The result is returned in a newly allocated matrix.
*/ public double[][] matmult( double matA[][], double matB[][] ) { // A is an MxN matrix int lenA_d1 = matA.length; int lenA_d2 = matA[0].length; // B is an NxM matrix int lenB_d1 = matB.length; int lenB_d2 = matB[0].length; double rslt[][] = null; if (lenA_d1 != lenB_d2 || lenA_d2 != lenB_d1) { System.out.println("matrix A must be MxN and matrix B must be NxM"); System.out.println("matrix A is " + lenA_d1 +" x " + lenA_d2 ); System.out.println("matrix B is " + lenB_d1 +" x " + lenB_d2 ); } else { int M = lenA_d1; int N = lenA_d2; rslt = new double[M][M]; for (int i = 0; i < M; i++) { for (int j = 0; j < M; j++) { double dotProd = 0.0; for (int k = 0; k < N; k++) { dotProd = dotProd + matA[i][k] * matB[k][j]; } rslt[i][j] = dotProd; } } } return rslt; } // matmult /**Matrix multiply
This is the linear algebra text book function to mltiply two matrices, where matA is MxN and matB is NxM, yielding an MxM result. The result is returned in argument matC, whose size must be MxM. */ public void matmult(double matC[][], double matA[][], double matB[][]) { // A is an MxN matrix int lenA_d1 = matA.length; int lenA_d2 = matA[0].length; // B is an NxM matrix int lenB_d1 = matB.length; int lenB_d2 = matB[0].length; // C should be an MxM matrix int lenC_d1 = matC.length; int lenC_d2 = matC[0].length; if (lenA_d1 != lenB_d2 || lenA_d2 != lenB_d1) { System.out.println("matrix A must be MxN and matrix B must be NxM"); System.out.println("matrix A is " + lenA_d1 +" x " + lenA_d2 ); System.out.println("matrix B is " + lenB_d1 +" x " + lenB_d2 ); } else if (lenA_d1 != lenC_d1 || lenA_d1 != lenC_d2) { System.out.println("matrix C should be " + lenA_d1 + " x " + lenA_d1 ); System.out.println("matrix C is " + lenC_d1 +" x " + lenC_d2 ); } else { int M = lenA_d1; int N = lenA_d2; for (int i = 0; i < M; i++) { for (int j = 0; j < M; j++) { double dotProd = 0.0; for (int k = 0; k < N; k++) { dotProd = dotProd + matA[i][k] * matB[k][j]; } matC[i][j] = dotProd; } } } } // matmult /** Transpose an MxN matrix into an NxM matrix */ double[][] transpose( double m[][] ) { int M = m.length; int N = m[0].length; double t[][] = new double[N][M]; for (int i = 0; i < M; i++) { for (int j = 0; j < N; j++) { t[j][i] = m[i][j]; } } return t; } // transpose /** Compare two vectors, a and b. Return true if the vectors are equal. */ boolean vecCompare( double a[], double b[] ) { boolean rslt = true; int N = a.length; if (b.length != N) { rslt = false; } else { for (int i = 0; i < N; i++) { if (a[i] != b[i]) { rslt = false; break; } } } return rslt; } // vecCompare /** Compare two matrices, a and b. Return true if they are equal, false otherwise. */ boolean matCompare( double a[][], double b[][] ) { boolean rslt = true; int M = a.length; int N = a[0].length; if (M != b.length || N != a[0].length) { rslt = false; } else { for (int i = 0; i < M; i++) { if (! vecCompare(a[i], b[i]) ) { rslt = false; break; } } } return rslt; } // matCompare } // linearOps matrix/linear/linearPr.java 100644 1040 1001 3567 7422400423 14456 0 ustar everyone package linear; /**
Print vectors and matrices
A vector can be printed in either vector form or as an ordered wavelet transform.
*/ public class linearPr { /** Round a double to four decimal places. */ private double round4( double d ) { final double multiplier = 10000.0; double rounded; rounded = Math.round(d * multiplier)/multiplier; return rounded; } // round4 /**Print a vector as an ordered wavelet transform. An ordered wavelet transform starts which a wavelet scaling value (sometimes called a smoothed value or an average). This is followed by bands of coefficients, where the length of the bands increases as a power of two (e.g., 20, 21, 22, ...)
@param vec[] result of a wavelet transform */ public void pr_ordered( double[] vec ) { if (vec != null) { int len = vec.length; if (len > 0) { System.out.println(round4(vec[0])); int num_in_freq = 1; int cnt = 0; for (int i = 1; i < len; i++) { System.out.print(round4(vec[i]) + " "); cnt++; if (cnt == num_in_freq) { System.out.println(); cnt = 0; num_in_freq = num_in_freq * 2; } } } } } // pr_ordered /** Print a vector */ public void prVec( double v[] ) { for (int i = 0; i < v.length; i++) { System.out.print( round4(v[i]) + " " ); } System.out.println(); } // prVec /** Print an MxN matrix */ public void prMat( double m[][] ) { if (m != null) { int M = m.length; int N = m[0].length; for (int i = 0; i < M; i++) { for (int j = 0; j < N; j++) { System.out.print( round4(m[i][j]) + " " ); } System.out.println(); } } } // prMat } matrix/matWavelet.java 100644 1040 1001 4771 7422403521 13541 0 ustar everyone import linear.*; /** Wavelet forward and inverse transform using matrices. */ class matWavelet { /** Move the odd elements to the second half of the N element region in the vec array. */ private void split( double[] vec, int N ) { int start = 1; int end = N - 1; while (start < end) { for (int i = start; i < end; i = i + 2) { double tmp = vec[i]; vec[i] = vec[i+1]; vec[i+1] = tmp; } start = start + 1; end = end - 1; } } /** Merge the odd elements from the second half of the N element region in the array with the even elements in the first half of the N element region. The result will be the combination of the odd and even elements in a region of length N. */ private void merge( double[] vec, int N ) { int half = N >> 1; int start = half-1; int end = half; while (start > 0) { for (int i = start; i < end; i = i + 2) { double tmp = vec[i]; vec[i] = vec[i+1]; vec[i+1] = tmp; } start = start - 1; end = end + 1; } } /**Calculate one wavelet transform step on the data set in "v" using a transform matrix m. The result is returned with the wavelet result (wavelet coefficients) in the upper half of the array and the scaling result (smoothed values) in the lower half.
The ordering is the matrix is scaling, wavelet, scaling, ... wavelet.
Note that the matrix m must be NxN
*/ public void transStep( double m[][], double v[] ) { int N = m.length; double scale, wave; linearOps ops = new linearOps(); for (int i = 0; i < N; i = i + 2) { scale = ops.dotprod(m[i], v ); wave = ops.dotprod(m[i+1], v ); v[i] = scale; v[i+1] = wave; } split( v, N ); } // transStep /** Calculate one inverse transform step on the input vector v using the inverse transform matrix m. The result is returned in the input vector v. */ public void invTransStep( double m[][], double v[] ) { int N = m.length; double v1, v2; linearOps ops = new linearOps(); merge(v, N); for (int i = 0; i < N; i = i + 2) { v1 = ops.dotprod( m[i], v ); v2 = ops.dotprod( m[i+1], v ); v[i] = v1; v[i+1] = v2; } } // invTransStep } // matWavelet