Continuous Wavelet Transforms, a Java Implementation

• Download source - 1.1 MB
• Download CwtRunnable.zip - 209 KB

sha256sum CwtSrc.zip: 1e3d72421b8d345a0b2f6f90cc209f8e6a8aa13f17c5bb6bc1a6fdbf5721729f

sha256sum CwtRunnable.zip: c3d872322f4b7803c22c7a1264dcd4493ae75eecc197c44e470d1d39e14223b6

Introduction

This is a companion to an earlier article I wrote on discrete wavelet transforms.  In this article I provide an application that uses continuous 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 present the underpinnings of the computations directly in Java methods.

Background

In the Fourier signal analysis, any signal may be understood as a linear combination of superimposed complex waveform components, each of which is stationary. The discrete Fourier transform of a signal reveals the frequency of the components that, resolution permitting, could be summed to re-create the original signal. But, because the components of the signal are stationary functions of frequency, all frequencies determined to exist in the signal are interpreted as existing throughout the entire duration of the signal.  What if the signal has localized structure that needs to be scrutinized, for example, temperature trends in climate histories, events in seismic records, crashes in an economy, etc.?

One solution is to iteratively perform a series of short term Fourier transforms over subsampled time windows within the signal. Another approach is to examine the convolution of the signal and a short term function - a wavelet - that can be iteratively scaled and translated. A CWT at scale s is the convolution of the signal with a translated version of the wavelet at scale s. Together the CWTs for multiple scales comprise a multi-resolution analysis (MRA).

 

One example of a wavelet: the 6^th derivative of the Gaussian distribution

My first article on discrete wavelet transforms (DWT) demonstrated this process for the discrete wavelet transform, which iteratively performs the transformation using successively halved wavelet scales. The continuous wavelet transform (CWT) allows for finer separations between adjacent scales. The term continuous is misleading as the analysis is performed on a discretely sampled signal, at discretely separated scales. For some applications the CWT reveals structure better than the DWT. For example, one of the Points of Interest in the DWT article was the revealed self-similarity of the Cantor staircase function. This is better demonstrated with the CWT.

Due to its relative infancy and opportunities for adaptation, wavelet analysis implementations suffer from non-standardized, sometimes-conflicting, application-specific treatments. The wavelets themselves have multiple synonyms. The parameters of the transform are expressed in different ways. The terms for characteristic quantities vary. Some treatments do not normalize scales while others use different (sometimes arbitrary) normalization approaches. The data presentation varies tremendously. To understand where the implementations in this article fit into this diaspora, the following reference is quite important and easy to get:

Christopher Torrence & Gilbert Compo, A Practical Guide to Wavelet Analysis, Bulletin of the American Meteorological Society, 20 October 1997.  

The development of my CWT implementations was originally based on study of their publication and algorithmic analysis of Fortran functions provided in their wavelet.f Fortran code. The above authors have gone beyond the scope of my article, particularly with respect to the statistical analysis required for meaningful interpretation of results in their demonstrated geophysical application.

The purpose of this article’s demonstration is to provide an entry point for programmers interested in using the vast power of the CWT. It provides an environment for exploring multi-resolution analysis of several functions, some of which are analytic and have verifiable revelations.  

Trying the runnable jar

To run the demo, get the CwtRunnable.zip download above, unzip it, find the CwtRunnable.jar file, and run it. On Windows, I simply double-click the CwtRunnable.jar file to run it. On Linux, I type java -jar <thePath>/CwtRunnable.jar into a terminal. Clearly, you will need to have Java installed, preferably the latest version with all updates.

Here is the intuitive opening form. First, try the defaults, then experiment.

The main entry point: Plot Generator

Clicking Make table will tabulate the function in a form. From here it can be plotted, transformed, saved, etc.

A function’s table form and its plot

Clicking on a plot will show the coordinates of the clicked point in the text box.

The table’s CWT menu will open a form for entering the transform wavelet and associated parameters. When you click OK, a multi-resolution table is created.

CWT parameter picker and MRA table

 

Now we can inspect the data graphically. First, select the Simultaneous Plots menu item. As the name suggests, it graphs the components of the multi-resolution analysis simultaneously on the same chart. Here is the simultaneous plot for the chirp example:

Simultaneous MRA plots for the chirp example

The Plots menu presents a picker table that you can use to selectively hide/show components. Zoom the x axis with [cntl] up and [cntl] down and the y axis with [alt] up and [alt] down. Again, clicking the chart shows the coordinates.

The Scalograms menu choice on the MRA table's menu provides the most common presentation given in the liturature.

Included in the unzipped demo are two comma-separated data files. From the plot generator form, click Open xy from csv and navigate to the TestSignal.csv file in your download folder. Plot the function from the data table form.

The signal was synthesized using the following components in a linear combination:

• From t = 0 to 0.255875,  dt = 0.000125
  • Component 1: sin(2pi500t)
  • Component 2: sin(2pi1000t)
  • Component 3: 0.1sin(2pi8t)
  • Component 4: a pulse with magnitude 5 located at t = 0.15 and 0.154

Obtain a CWT and select Scalogram > Modulus from the MRA table’s menu. The signal’s component specification tells us that we should have primary dark bands near f = 8, f = 500, and f = 1000.  We should also have two impulses at 0.150 and 0.154.

As you try the examples, some scalograms will require contrast adjustment, a feature controlled by the mouse wheel.

Choosing Reconstruction from the MRA table form's menu provides a reality check for the applicability of your chosen wavelet and parameter choices for effectively analyzing your signal. For maximal reproduction, use a low value for the minimum scale and a high value for the total number of scales on the parameter input form. A low relative residual indicates a good reconstruction. Although lowering the smallest scale and increasing the number of scales in the transform improves the reconstruction, it does not necessarily increase the usability of the MRA data itself. Here we must consider the cone of influence.

The CWT implementation involves iteratively convolving scaled wavelets with the signal using the convolution theorem. The convolution is the inverse Fourier transform of the point-wise product of the Fourier transforms of the two components. Finite, discrete samplings of signals have abrupt boundaries at the start and end of the sequence and convolutions that overlap the boundary are contaminated with errors. The effect can be visualized with a function that has a pulse at the beginning and end of the signal and zero elsewhere. Here is the MRA scalogram of such a function.

The grey area proportionally shows the regions of the MRA that are influenced by the error (within the cone of influence).  Transforms at wider wavelet scales (towards the top) lie increasing within the grey area. There are various ways to graphically compensate for the error (e.g. cosine damping) but they require artificial distortions/extrapolations of the signal. The scalogram of the Devil’s staircase function at the beginning of this article was purposely cosine damped to produce the graphic. When possible one can collect the signal over a longer period of time and selectively ignore features early and late in the MRA. There are other ways to manage the error, such as symmetric extension, but the applicability is application specific.

Unlike synthetic functions, signals arising from natural systems are noisy. One must statistically determine, to a pre-defined confidence level, the significance of the MRA observations.The question asked is: are we reasonably confident that an observed feature is not actually a noise excursion? Answering the question requires an understanding of how noise distributes in the system being studied and a knowledge of related statistics. This is application-specific and, because this demo focuses on exact synthetic signals,  significance testing is not addressed.

Using the code

A license is included in the header for each class: Copyright Mark Bishop, GNU v3.  Credits to others are cited. The author makes no warranty for the accuracy, completeness, safety, or usefulness of any information provided and does not represent that its use would not infringe privately owned right.

The CwtSrc.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 CwtDemo. 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 to use the example function expressions. It is located in the jel-2.0.1 folder. I submitted an article on using the JEL expressions library at www.codeproject.com/Articles/737497/Java-Function-Compiler. Jel expressions are a part of the demo but are not required for using the CWT specific classes.

The computational classes are the only classes required for obtaining CWT multi-resolution analyses. The remaining classes and the JEL expressions library are ancillary features of the demo that provide a graphical interface with tables, dialogs, plots, etc.

The following classes are included in the project:

• Computational classes
  • ComplexCalc.java (Static methods for complex calculations)
  • ComplexNumber.java (Object for containing a single complex number)
  • CWT.java (The primary subject of the article)
  • Fft.java (Forward and reverse Fourier transform methods)
  • Gamma.java (Computes the gamma function)
  • Hermite.java (Methods for computing Hermite polynomials)
  • MatrixOps.java (Used here for array management)

• Interactive GUI classes
  • MainForm.java (The main plot generator form)
  • MraForm.java (The MRA table form)
  • XYForm.java (The function's table form, t and f(t))
  • PlotFrame.java (The function and MRA graphs)
  • PlotVisiblePicker.java (A visibility selector for simultaneous plots)
  • ParameterInputCwt.java (Dialog form for entering transform paramters)
  • Scalogram.java (Scalogram presentation)
• Other utilities
  • Examples.java (A few example functions of interest)
  • FileOps.java
  • ImagePanel.java
  • Signal.java (a helper for managing xy data)
  • StringUtils.java
  • ConsoleTest.java (Limited reproduction testing routine for each wavelet)
  • VariableProvider.java (required by jel)

The class of primary importance is CWT.java. Within this class, a single method is responsible for returning the multi-resolution data: cWT().

/**
     * Continuous wavelet transform for a real, discretely sampled function y =
     * f(t) of length n where n is an integer power of 2
     *
     * @param y
     *            The time-sampled signal to transform.
     * @param dt
     *            The sampling increment
     * @param mother
     *            The wavelet to use in the transform (Morlet, DOG, Paul)
     * @param param
     *            wavelet parameter (Morlet param = wave number, Paul param =
     *            order, DOG param = m::m-th derivative of Gaussian)
     * @param s0
     *            The smallest desired scale in the transform (try s0=dt for
     *            Morlet; s0=dt/4 for Paul)
     * @param dj
     *            The desired spacing between successive scales (try 0.25)
     * @param jtot
     *            The desired number of scales to be computed
     * @return An ArrayList[Object] wspcmf[6] where: xspcmf[0] is a
     *         Complex[n][jtot] multi-resolution matrix, wspcmf[1] is a
     *         double[jtot] vector of scales, wspcmf[2] is a double[jtot] vector
     *         of periods, wspcmf[3] is a double[n] cone-of-influence vector,
     *         wspcmf[4] is the signal mean, and wspcmf[5] is the Fourier
     *         factor.
**/		
public static ArrayList<Object> cWT(double[] y, double dt, Wavelet mother, double param, double s0, double dj, int jtot)

• Input parameters:
  • y[]: The signal for transformation
  • mother: The wavelet (Morlet, derivative-of-Gaussina (DOG), or Paul)
  • param: The wavelet parameter (Morlet: wave number, Paul: order, DOG: m where we use mth derivative of the Gaussian distribution)
  • s0: the smallest scale to be included in the MRA
  • dj: The spacing between scales
  • jtot: The total number of scales requested

• Output: An ArrayList<Object> named wspcmf with the following contents
  • xspcmf[0] is a ComplexNumber[n][jtot] multi-resolution matrix
  • wspcmf[1] is a double[jtot] vector of scales
  • wspcmf[2] is a double[jtot] vector of periods
  • wspcmf[3] is a double[n] cone-of-influence vector (unused in demo)
  • wspcmf[4] is the signal mean
  • wspcmf[5] is the Fourier factor

The method first creates a copy of y, named yNoDC, with the signal mean removed. Next, its Fourier transform is obtained as yfft. After creating a vector of wave numbers, the main loop calls the method getDaughter(), which returns the Fourier transform of the appropriate wavelet function at the scale of iteration (see Torrence & Compo, Table 1). The loop then performs the convolution at the scale by obtaining the inverse Fourier transform of the product of the daughter and yfft.

The CWT at each scale can be tabulated or co-plotted vs. the set of scales, the periods, or the reciprocals of the periods (frequencies). In the demo the reciprocal periods are used. The Fourier factor and the period vector are included in the return object for convenience and could alternately be computed as necessary in a consuming method. The cone-of-influence may be processed to envision the areas where signal boundary effects contaminate the MRA.

The cwtReconstruct() method of CWT.java reconstructs the original signal from the real part of the MRA matrix and the vector of scales. It requires the original signal’s mean value that was subtracted from the signal during the transform as well as the sampling increment, dt. The other input parameters must match the values used in the original transform. The cwtReconstruct() method calls the getEmpiricalFactors() method, which determines the  reconstruction factors for the wavelet with a delta function calibration (see Torrence & Compo). Have a look at my article, entitled  Continuous Wavelet Transform Reconstruction Factors for Selected Wavelets, located at: mark-bishop.net/signals/CWTReconstructionFactors.pdf. Clearly, the reconstruction factors could have been hard-coded rather than computed each time a reconstruction is performed.

The cWTReconstruct.java class also includes a convenience method, getSelectedParamChoices(), which provides some of the possible wavelet parameters for the three wavelets in the demo. The DOG 5 choice is included as a brain teaser.  Finally, the class has an enumeration, Wavelet, containing the three wavelet choices.

A few references

Christopher Torrence & Gilbert Compo, A Practical Guide to Wavelet Analysis, Bulletin of the American Meteorological Society, 20 October 1997.

Marie Farge, Wavelet Transforms and Their Applications to Turbulence, Annu. Rev. Fluid Mechanics, 24:395-457, 1992

S. D. Meyers, Kelly, Obrien; An Introduction to Wavelet Analysis in Oceanography and Meteorology: With Application to the Dispersion of Yanai Waves, Mon. Wea. Rev., 121, 2858-2866, 1993.

Stéphane G. Mallat, A Wavelet Tour of Signal Processing, The Sparse Way, 3rd ed, Academic Press 2008.

E. Oran Brigham, The Fast Fourier Transform, Prentice-Hall, 1974.

Paul Addison, The Illustrated Wavelet Transform Handbook, Taylor & Francis, 2002.

Jaideva Goswami, Andrew Chan, Fundamentals of Wavelets, Wiley, 1999.

http://users.rowan.edu/~polikar/WAVELETS/WTtutorial.html

 

History

 

Keep a running update of any changes or improvements you've made here.