DwtSrc.zip sha256 sum: 665bfbd806214149bd569571255e1cc15ffd660124cede54caf22bf7964db86c
DwtRunnable.zip sha256 sum: 95030fe0f023bcefc9874c5580b42dddb6acfe0f14bc9caf3d6c635f5b844bd4
Introduction
Wavelet transforms are essential tools for making sense of the world (AKA signal processing). The term signal refers to any information stream that varies as a function of some range variable(s). Most often, in signal processing literature, the range variable is time, but it could be essentially anything: position, applied-magnetic field, wavelength, or even a dimensionless sequence. In applications, the function can be a recorded sound, an image, a transducer response, an EKG or MRI, spectrometry data, climate records, chaos, etc. The list is endless and there are countless creative opportunities for developers to explore novel approaches to data analysis using wavelets. One very active dicipline is machine-based pattern recognition.
In this article I provide an application that uses discrete wavelet transforms to explore one dimensional signals. The approach can be further developed to transform signals with higher dimensions, like images. There are other ways to add wavelet functionality to Java, such as employing Matlab/Scilab-Wavelab wrappers or open source libraries. For illustration, I chose to present the underpinnings of the computations directly in Java methods.
The computational classes in this project are translations of C# code I wrote a couple of years ago. Reversing the process should be fairly easy.
Background
You may have arrived at this article after searching Google for something like: compute wavelet transform. Some of the hits provide important definitions and mathematical background with varying degrees of complexity. Others answer the question of how to compute a transform by simply stating something like: open Matlab/Scilab/etc and type...
Those answers are valid but developers want to see, or work out, an implementation of the algorithms in their language of choice. Below are mine. There are more elegant and efficient approaches but I chose to literally express the algorithms in methods that run on a single CPU thread.
One challenge in developing this demo was effective visualization of multi-resolution data. Printed presentations require the use of stacked plots, which are difficult to evaluate interactively on a computer monitor. In this case I chose to simultaneously plot the scales on the same graph, allowing the user to show/hide individual scales from a menu item. Also, I structured the data so that the superimposition of all scale data would perfectly reconstruct the signal.
Trying the runnable jar
To run the demo, download the DwtRunnable.zip download above, unzip it, find the DwtRunnable.jar file, and run it. On Windows, I simply double-click the DwtRunnable.jar file to run it. On Linux, I type java -jar <thePath>/DwtRunnable.jar into a terminal. Clearly, you will need to have Java installed, preferably the latest version with all updates. There are examples you can run from the main GUI. Click Plot to see the graph and the signal's data table. From the data table window you can perform transforms or multi-resolution analysis. You can load the menu examples or you can type a function into the expression box. Also, if you click Empty Table you will be presented with a blank table where you can open a csv file for plotting and analysis.
Here are some screen shots using the Gauss modified mixed sin example.
Here is the initial GUI, its plot, and its data table:
Here is the transform input dialog, the discrete wavelet transform, and its inverse (reconstruction):
Finally, here's the multi-resolution analysis and its table:
The columns in the multi-resolution table will sum to the original signal. The individual scales in the MRA plot can be shown/hidden from the menu.
Using the code
The DwtSrc.zip download is a zipped Eclipse project folder. If you use Eclipse, extract it into the workspace of your choosing and, from the IDE, create a project exactly named DwtDemo. This should populate the project and enable a test run. With other IDE's or if you choose another way for creating your Eclipse project, you will need to link to the jel.jar library in your build path. It is located in the jel-2.0.1 folder. There is a previous article on using the JEL expressions library at http://www.codeproject.com/Articles/737497/Java-Function-Compiler. It is part of the demo but is not required for using the discrete wavelet transform specific classes.
The following 16 classes are included:
- Colors.java (a few colors that plot multi-resolution graphs nicely)
- Wavelet.java (Haar, Daubechies, Symmlet, Coiflet, and others)
- JFrame based interactive GUIs
- MainForm.java (entry point for demo)
- MatrixForm.java (data tables and related menus)
- ParameterInput.java (user selection of transform parameters)
- PlotFrame.java (plots with menus)
- DWT.java (primary subject of article)
- OrthogonalFilters.java (also primary)
- MatrixOps.java (used by DWT.java)
- VisualizeMultiResolution.java (application specific helper)
- Examples.java (a few special functions, some non-analytic)
- FileOps.java
- ImagePanel.java
- Signal.java (application specific helper for managing xy data)
- StringUtils.java
- VariableProvider.java (required by JEL)
You will find a comment near the beginning of each class that states the class's responsibility. I tried to properly comment the methods where needed. Error handling is minimal. There are a few important things to discuss about DWT.java and OrthogonalFilters.java.
DWT.java includes the method transform():
public static double[] transform(double[] signal, Wavelet wavelet,
int order, int L, Direction direction) throws Exception {
if (direction == Direction.forward) {
return forwardDwt(signal, wavelet, order, L);
} else {
return inverseDwt(signal, wavelet, order, L);
}
}
- The signal is the sequence you want to transform, which must have a length that is an even power of 2. (Another method is provided for padding if needed.)
- Wavelet can be Haar, Daubechies, etc.
- Order is a wavelet specific choice that must agree with one of the choices defined for each wavelet type provided in OrthogonalFilters.java. The nomenclature and form of these filter designations varies considerably in the literature.
- L is the coarsest scale that you want to include in your transform.
- Direction is forward for transforms and reverse for inverse transforms (defined in the class).
Depending on direction, it will call either forwardDwt() or inverseDwt().
For the forward transform, I prefer to iteratively form, for each scale, a quadrature mirror filter (QMF) as a matrix using the appropriate high-pass and low-pass filters provided in OrthogonalFilters.java. The mirror matrices are created by makeQMFMatrix() at each scale and a cascading matrix-vector multiplication implements the transform. The return result, dWT, is initially set to the input signal. Each time the scale is halved, a vector (subResult) containing the first half of dWT is multiplied by the appropriately sized mirror matrix. Subsequently, the product overwrites the first n/(2i) values of dWT, where i is the current iteration in the for loop.
public static double[] forwardDwt(double[] signal, Wavelet wavelet,
int order, int L) throws Exception {
int n = signal.length;
if (!isValidChoices(wavelet, order, L, n)) {
throw new Exception(
"Invalid wavelet /order/scale/signal-length combination.");
}
double[] dWT = MatrixOps.deepCopy(signal);
int log2n = (int) (Math.log(n) / Math.log(2));
int iterations = log2n - L;
int subLength = n;
double[] H = OrthogonalFilters.getLowPass(wavelet, order);
double[] G = OrthogonalFilters.getHighPass(H);
for (int i = 0; i < iterations; i++) {
subLength = n / (int) (Math.pow(2, i));
double[][] QMF = makeQMFMatrix(subLength, H, G);
double[] subResult = new double[subLength];
subResult = subCopy(dWT, subResult, subLength);
double[] temp = MatrixOps.multiply(QMF, subResult);
dWT = subCopy(temp, dWT, subLength);
}
return dWT;
}
The inverseDwt() method can be applied to any signal, but usually it is applied to a previously obtained forward transform for the purpose of reconstruction, often after some modification, like de-noising, is performed on the transform. If you are reconstructing, use the identical wavelet, order, and coarsest scale (L) that you used in the original transform. The method reverses the actions of a forward transform in a manor similar to the forwardDwt(), except that the transpose of the QMF matrix is used as a multiplier within a for (i = L + 1; i <= log2n, i++) loop structure.
public static double[] inverseDwt(double[] signal, Wavelet wavelet,
int order, int L) throws Exception {
int n = signal.length;
if (!isValidChoices(wavelet, order, L, n)) {
throw new Exception(
"Invalid wavelet /order/scale/signal-length combination.");
}
int log2n = (int) (Math.log(n) / Math.log(2));
int subLength;
double[] preserveCopy = new double[signal.length];
preserveCopy = subCopy(signal, preserveCopy, signal.length);
double[] H = OrthogonalFilters.getLowPass(wavelet, order);
double[] G = OrthogonalFilters.getHighPass(H);
for (int i = L + 1; i <= log2n; i++) {
subLength = (int) (Math.pow(2, i));
double[][] QMF = makeQMFMatrix(subLength, H, G);
QMF = MatrixOps.transpose(QMF);
double[] subResult = new double[subLength];
subCopy(signal, subResult, subLength);
subResult = MatrixOps.multiply(QMF, subResult);
signal = subCopy(subResult, signal, subLength);
}
double[] iDWT = new double[n];
iDWT = subCopy(signal, iDWT, n);
signal = preserveCopy;
return iDWT;
}
There are many places to schematically review the underlying algorithms implemented above including: http://users.rowan.edu/~polikar/WAVELETS/WTpart4.html figure 4.1. The QMF matrix implementation used here is more difficult to find, so convince yourself that it is functionally correct (some tests are provided below). The demo's implementation is somewhat inefficient because, as implemented, it explicitly processes a sparse matrix as if it were a dense matrix. The conceptual model, filtering using a matrix-vector product, is very intuitive and could be re-implemented using a sparse multiplication strategy, perhaps on parallel GPU threads or on a GPU.
Multi-resolution analysis
The method that performs multi-resolution analysis, mRA(), takes the same set of arguments as the transform() method and returns an ArrayList<Object> with the following specification
- mRA.get(0) is an ArrayList<double[]> holding multi-resolution scale data from the finest scale to the coarsest scale, plus the final approximation scale that contains the curve corresponding to scales greater than L.
- mRA.get(1) is an int[] containing values j, such that each scale in the MRA, from finest to coarsest, may be determined later with the formula: scale = 2-j
- The final double[] listed in mRA.get(0) and the final j provided by mRA.get(1) correspond to the approximation curve, which carries frequencies larger than 2-L
Below is the mRA() method. Initially the signal is forward-transformed into dwt. The first loop creates the first n-1 items in the ArrayList<double[]> named mRA by zero-interpolating successive diadic scales in dwt, then inverse-transforming the intermediate result. The second loop adds the final approximation curve for frequencies larger than 2-L . The third loop adds the j values to an int[] named scalesUsed, designating the final approximation curve with j = 0. Finally, mRA and scalesUsed are added as Objects to the return value, named result.
public static ArrayList mRA(double[] signal, Wavelet wavelet,
int order, int L) throws Exception {
ArrayList result = new ArrayList();
int n = signal.length;
if (!isValidChoices(wavelet, order, L, n)) {
throw new Exception(
"Invalid wavelet /order/scale/signal-length combination.");
}
int J = (int) (Math.log(n) / Math.log(2));
double[] dwt = forwardDwt(signal, wavelet, order, L);
ArrayList mRA = new ArrayList();
for (int j = (J - 1); j >= L; j--) {
double[] w = new double[n];
int[] dyad = dyad(j);
for (int k = dyad[0]; k <= dyad[dyad.length - 1]; k++) {
w[k - 1] = dwt[k - 1];
}
mRA.add(inverseDwt(w, wavelet, order, L));
}
double[] w = new double[n];
int limit = (int) Math.pow(2, L);
for (int i = 0; i < limit; i++) {
w[i] = dwt[i];
}
mRA.add(inverseDwt(w, wavelet, order, L));
int[] scalesUsed = new int[mRA.size()];
int scaleCounter = 0;
for (int j = (J - 1); j >= L; j--) {
int[] dyad = DWT.dyad(j);
scalesUsed[scaleCounter] = (int) (Math.log(dyad.length) / Math
.log(2));
scaleCounter++;
}
scalesUsed[scaleCounter] = 0;
result.add(mRA);
result.add(scalesUsed);
return result;
}
Code testing
I tested my results in several ways.
First, I transformed several functions using various wavelets and wavelet parameters, inspected the result to be certain it conformed with expectations, then inverse transformed the result and verified perfect reconstruction. In the DemoSrc download directory there is a spreadsheet named TestSignal.xls that presents a 3-frequency, mixed sinusoid function with momentary transient noise between x = 0.0525 and x = 0.0625, and two spikes at 0.154 and 0.155. This signal was used to create a csv named TestSignalXY.csv.
There is a test method in DWT.java named relativeErrors() that enables automated reconstruction testing for all orthogonal filter combinations provided in OrthogonalFilters.java. Running this method with TestSignalXY.csv as a data file shows relative residuals for the difference between the original signal and the reconstructed signal for each wavelet/order combination. The values range between ~10-12 to ~10-16 for all wavelets except the Battle wavelet, which is ~10-4. The coefficients for the Battle wavelet are only provided to 9 significant figures and reconstruction is expected to be diminished. Using the application, the columns of the MRA matrix add to perfectly reproduce the signal. Multi-resolution plots for selected scales are shown below in stacked format. The finest scale shows the spikes, the mid-scale shows the transient, and the remaining plot shows the three sinusoids.
I also analyzed several test signals using Scilab/Wavelab and verified that the demo's results for forward and reverse transforms agree. A spread sheet is included in the download that presents one case.
Points of Interest
Try the Dropout in signal example from the menu. Obtain a multi-resolution plot from the transform menu of function's data table. From the MRA's plot menu, hide all scales except the 2-9 scale. This demonstrates how remarkably well the wavelet transform detects the dropout flaw.
The point where the flaw occurs corresponds to modulus maximum at the scale, perhaps indicating a singularity in the signal. In real-world examples, maxima may represent boundaries in an image or a crack in a machine part, revealed by a vibration profile.
A random sequence signal would (most often) be a singular function, displaying modulus maxima throughout each scale.
The Cantor function is a singular function sometimes called the devil's staircase. The scales of its transform show its self-similarity nicely.
A few references
Stéphane G. Mallat, A Wavelet Tour of Signal Processing, The Sparse Way, 3rd ed, Academic Press 2008
http://statweb.stanford.edu/~wavelab/Wavelab_850/Contact.html.
Paul Addison, The Illustrated Wavelet Transform Handbook, Taylor & Francis, 2002.
Jaideva Goswami, Andrew Chan, Fundamentals of Wavelets, Wiley, 1999.
Martin Vetterli, Wavelets and Filter Banks: Theory and Design, IEEE Trans on Sig Proc, Vol 40 No 9, 1992.
http://waveletsandsubbandcoding.org/Repository/VetterliKovacevic95_Manuscript.pdf
http://users.rowan.edu/~polikar/WAVELETS/WTtutorial.html
JEL: http://www.codeproject.com/Articles/737497/Java-Function-Compiler
History
None yet.
