Spline Wavelet Manual

Susumu Sakakibara
Iwaki Meisei University

May 1994

This Notebook explains how to use the Mathematica Package,
SplineWavelet.m, developed by the author.

Theory

A scaling function phi[x] is a square integrable function which satisfied the two-scale relation

phi[x] = Sum[p[[l]] phi[2 x - l], {l, ...}]

where p is a square summable sequence. Usually phi[x] is normalized such that

Integrate[phi[x], {x, -Infinity, Infinity}] = 1

Let Vj be the space of all functions that can be written as

f[j, x] = Sum[c[[l]] phi[2^j x - l], {l, ...}]

If f[x] belongs to Vj, then f[2 x] belongs to Vj+1. Thus the spaces have the hierarchical structure that Vj-1 is contained in Vj, which is contained in Vj+1, and so on. This structure is called the multiresolution analysis (MRA) generated by the scaling function phi[x].

Let Wj be the complementary subspace of Vj, and let psi[x] be a function such that psi[x - l], where l are integers, form a basis of W0, with

Integrate[psi[x], {x, -Infinity, Infinity}] = 0

The function psi[x] is called the wavelet. Any function g[j, x] belonging to Wj may be written as

g[j, x] = Sum[d[[l]] phi[2^j x - l], {l, ...}]

Due to the structure Vj+1 = Vj + Wj, any function f[j, x] belonging to Vj can be uniquely written as

f[j, x] = f[j - 1, x] + g[j - 1, x]

If we find f[0, x] that fits the data sequence, then we can apply this decomposition, repeatedly if needed, to get f[-1, x], f[-2, x], and so on. Since psi[x] has an oscillatory behavior, this decomposition can be used to smooth out the data, thereby eliminating noise.

Chui and Wang (see C. K. Chui, An Introduction to Wavelets, Academic Press, 1992) have constructed scaling functions and wavelets with compact supports, based on the cardinal B-splines. The package SplineWavelet.m defines the Mathematica functions that can be used to data smoothing described above. They are implemented by the algorithms using the cubic cardinal B-splines.

Initialization

Off[General::spell];
Off[General::spell1];
Needs["SplineWavelet`"]

Here it is assumed that the package SplineWavelet.m is located in the same level as the Mathematica kernel. If it is not the case, the loading path must be set appropriately before loading the package SplineWavelet.m.

The structure of the data

We shall deal with a set of data in the form of a nested list.

data = {-2, {1.2, 3.2, 2.8, 1.6, 1.2, 3.6, 0.4}};

The second element {1.2, 3.2, 2.8, 1.6, 1.2, 3.6, 0.4} is the data sequence,
fk, k = -2, -1, ..., 4, and the first element -2 denotes the first index of the data sequence. The values of fk that are not specified are assumed to vanish. The set of data can be plotted by the PlotData function.

pd = PlotData[data];

The scaling function

The function FineDisplayScaling plots the linear combination

f[j, x] = Sum[c[[l]] phi[2^j x - l], {l, ...}]

Thus by taking c0 = 1, with all other cl vanishing, we get a plot of phi[x] itself. The sequence c is given as {0, {1}}, as explained in the previous section.

FineDisplayScaling[{0, {1}}];

The function FineDisplayScaling implements the interpolatory graphical display algorithm and uses ListPlot to generate the above plot, and hence it is faster than Plot which uses adaptive sampling. To see how the algorithm works, let us use the DisplayScaling function instead.

DisplayScaling[{0, {1}}];

The plot is a ListPlot with the data points at the integer knots. These data points are obtained by the function ScalingSample.

ScalingSample[{0, {1}}]

The data points at the integer knots are not enough to show the function in detail. The interpolatory graphical display algorithm allows us to get the function values at dyadic points x = 2^-k n. This is done by specifying the second argument in the ScalingSample function. In general ScalingSample[data, {j, k}] generate the set of pairs {x, f[j, x]}. The parameter k specifies how many times the algorithm is iterated, so that k mus be a nonnegative integer. In the shorter form ScalingSample[data, j], k is set equal to 0. Thus the following command generates the list with x values at half-integer knots.

ScalingSample[{0, {1}}, {0, 1}]

The function FineDisplayScaling generates plots using this function with k = 4, using ListPlot with PlotJoined -> True by default. When the number of data points is very large, DisplayScaling would be sufficient to plot the function. In both of FineDisplayScaling and DisplayScaling, one can set PlotJoined -> False to override the dafault setting True.

Note that ScalingSample generates {x, y} pairs. One might just want a list of y values, e.g., to feed in a ParametricPlot, or to generate sound with ListPlay. For such purposes, one can use ScalingInterpolate instead. Thus

ScalingInterpolate[{0, {1}}]

generates the list of y values corresponding to ScalingSample[{0, {1}}]. If we spacify an optional interpolation level, we get, for example,

ScalingInterpolate[{0, {1}}, {1}]

Wavelet

The function FineDisplayWavelet plots the linear combination

g[j, x] = Sum[d[[l]] psi[2^j x - l], {l, ...}]

Thus by taking d0 = 1, with all other dl vanishing, we get a plot of psi[x] itself. The sequence d is given as {0, {1}}.

FineDisplayWavelet[{0, {1}}, PlotRange -> All];

This function also uses the interpolatory graphical display algorithm. In fact, the "original data" to produce this plot are the data at half integer knots, which may be seen by the function DisplayWavelet.

DisplayWavelet[{0, {1}}, PlotRange -> All];

These data are obtained by the function WaveletSample.

WaveletSample[{0, {1}}]

As the data points at the half-integer knots do not show the function in detail, we will use the interpolatory graphical display algorithm to get the function values at dyadic points x = 2^-k n. This is done by specifying the second argument in WaveletSample. In general, WaveletSample[data, {j, k}] generate the set of pairs {x, g[j, x]}. The parameter k - 1 specifies how many times the algorithm is iterated, so that k must be a positive integer. In the shorter form ScalingSample[data, j], k is set equal to 1. Thus the following command generates the list with x values at quarter-integer knots.

WaveletSample[{0, {1}}, {0, 2}]

The function FineDisplayWavelet generates plots using this function with k = 4, using ListPlot with PlotJoined -> True by default. When the number of data points is very large, DisplayWavelet would be sufficient to plot the function. In both of FineDisplayScaling and DisplayScaling, one can set PlotJoined -> False to override the dafault setting True.

The function WaveletInterpolate is similar to WaveletSample, but only produces y values.

WaveletInterpolate[{0, {1}}]

WaveletInterpolate[{0, {1}}, {2}]

If we use the flat data

flatData = {0, {1, 1, 1, 1, 1}};

for the scaling function and wavelet, we can see the basic properties of the scaling function and wavelet.

FineDisplayScaling[flatData];

FineDisplayWavelet[flatData];

The scaling function interpolates the data, while the wavelet has the oscillatory behavior.

Interpolation

We now want to find the scaling function f[0, x] that fits the data given at the beginning. We need to find the list of coefficients c[0], which can be found by the function FindInterpolation.

c[0] = FindInterpolation[data]

We can display the corresponding scaling function f[0, x] together with the original data.

scalingPlot =
	FineDisplayScaling[c[0], PlotRange -> All,
		DisplayFunction -> Identity];
Show[scalingPlot, pd, DisplayFunction -> $DisplayFunction];

Since the space V0 is contained in V1, f[0, x] may also be expressed in terms of the basis of V1.

f[0, x] = Sum[c[0][[l]] phi[x - l], {l, ...}]
		= Sum[c[1][[l]] phi[2 x - l], {l, ...}]

The function FinerScalingCoefficient generate the sequence c[1] from the sequence c[0].

c[1] = FinerScalingCoefficient[c[0]]

Hence the following command will generate the same plot, generated by FineDisplayScaling.

ListPlot[ScalingSample[c[1], {1, 3}], PlotJoined -> True];

Decomposition

To see how the decompostion works, we need a lot of data points. So, we will generate a noisy data as follows.

newData =
	{0, Table[N[Sin[(1 + Random[ ]/8) x]],
		{x, 0, 6 Pi, Pi/15}]};

newDp = PlotData[newData];

The idea here is to produce a noise, that are not a mixture of simple harmonics.

Using FindInterpolation, we get the sequence c[0].

c[0] = FindInterpolation[newData];

The fitting curve is then given by

dataPlot =
	DisplayScaling[c[0],
		PlotStyle -> {GrayLevel[0.4], Thickness[0.007]},
		DisplayFunction -> Identity];
Show[dataPlot, newDp, DisplayFunction -> $DisplayFunction];

This curve can be decomposed into a smoother part f[-1, x] and noisy part g[-1, x]. For the scaling function f[-1, x], the cofficient sequence is obtained as

c[-1] = DecomposeToScaling[c[0]];

We can plot f[-1, x] now.

DisplayScaling[c[-1], -1];

Here the second argument -1 in DesplayScaling specifies the resolution level; without this argument, we would have got the x axis half scaled.

The wavelet component can be obtained by the following way.

d[-1] = DecomposeToWavelet[c[0]];
juggy1 = DisplayWavelet[d[-1], -1, PlotRange -> All];

We can iterate the decomposition further.

c[-2] = DecomposeToScaling[c[-1]];
smoother = DisplayScaling[c[-2], -2, PlotRange -> All];

d[-2] = DecomposeToWavelet[c[-1]];
juggy2 = DisplayWavelet[d[-2], -2, PlotRange -> All];

Let us compare the smoother curve with the original data.

Show[dataPlot, smoother,
	DisplayFunction -> $DisplayFunction];

The sum of the three components is equal to the original function.

f[0, x] = f[-2, x] + g[-1, x] + g[-2, x]

which would correspond, in the graphical form, to

Show[GraphicsArray[{{smoother, juggy1, juggy2}}]];

Data smoothing

Note: To actually run the session of this section requires about 10MB of memery for Mathematica kernal.

We now investigate data smoothing with a realistic example. Let us read in the data of 6320 numbers.

rData = {0, ReadList["exp136b.03", Number]};

rDataP = PlotData[rData, PlotStyle -> {PointSize[0.005]},
	AspectRatio -> 1/4];

The coefficients for the interpolation is

c[0] = FindInterpolation[rData];

We will perform the decomposition twice.

c[-1] = DecomposeToScaling[c[0]];
sc1 = DisplayScaling[c[-1], -1, AspectRatio -> 1/4];

d[-1] = DecomposeToWavelet[c[0]];
wl1 = DisplayWavelet[d[-1], -1, PlotRange -> {-300, 300},
	AspectRatio -> 1/4];

If we compare the plots, we see that the noise is removed from the data.

Show[rDataP, PlotRange -> {{0, 500}, {-60, 60}}];
Show[sc1, PlotRange -> {{0, 500}, {-60, 60}}];
Show[wl1, PlotRange -> {{0, 500}, {-60, 60}}];

Show[rDataP, PlotRange -> {{2001, 2500}, {-350, 350}}];
Show[sc1, PlotRange -> {{2001, 2500}, {-350, 350}}];
Show[wl1, PlotRange -> {{2001, 2500}, {-350, 350}}];

It is remarkable that human ear is very sensitive to a periodic signal. In fact one can detect noise contained in the original data by generating sound.

ListPlay[rData[[2]], SampleRate -> 4000];

To generate sound for smoother function f[-1, x], we use ScalingInterpolate.

soundSc = ScalingInterpolate[c[-1], {1}];

soundWl = WaveletInterpolate[d[-1], {1}];

ListPlay[soundSc, SampleRate -> 4000];

ListPlay[soundWl, PlayRange -> {-300, 300},
	SampleRate -> 4000];

The data smoothing carried out here is not the one that would be obtained by a simple application of the Fourier method. To see this we plot the Fourier spectra of the original as well as the smoothed data.

ListPlot[Abs[Fourier[rData[[2]] ] ] ];

ListPlot[Abs[Fourier[soundSc ] ] ];

References

C. K. Chui, An Introduction to Wavelets, Academic Press, 1992.
S. Sakakibara, A Practice of Data Smoothing by Cardinal B-Spline Wavelets, 
     in Wavelets: Theory, Algorithms, and Applications, ed. by C. K. Chui, L. Montefusco,
     and L. Puccio, Academic Press, 1994, to appear.