How to define a new wavelet type?

Question

I have a problem in understanding how to create my own wavelets in Mathematica's wavelet framework.

There is a tutorial.

I have read the tutorial, but I don’t understand what exactly is needed to define my own wavelet type. I tried to implement a very simple wavelet (Haar wavelet) in the framework discussed in the tutorial to get a basic understanding, but it did not work out.

Does anyone has experience with this framework or directly sees how or what I need to do? I want to use this framework because it is the only way of using internal Mathematica functions with your own wavelets.

The tutorial says the following:

FranklinWavelet[lim_Integer]["WaveletQ"] := True
FranklinWavelet[lim_]["OrthogonalQ"] := True

a[k_] := Sqrt[3] HypergeometricPFQRegularized[{1/2, 1/2, 1}, {1 - k, 1 + k}, -2]

 Z[k_, n_] := 1/24 a[2 k - n - 2] + 1/4 a[2 k - n - 1] + 5/12 a[2 k - n] + 1/4 a[2 k - n + 1] + 1/24 a[2 k - n + 2]

FranklinWavelet[lim_Integer]["PrimalLowpass", prec_ : MachinePrecision] := 
 ParallelTable[{n, N[Total[Table[a[k] Z[k, n], {k, -15, 15}]], prec]}, {n,-lim, lim}]

wave = FranklinWavelet[14];

WaveletFilterCoefficients[wave]

{{-14, -9.77137*10^-6}, {-13, 0.0000212211}, {-12, 0.0000391428}, {-11, -0.0000867523}, {-10, -0.000158601}, {-9, 0.000361781}, {-8, 0.000652922}, {-7, -0.00155701}, {-6, -0.00274588}, {-5, 0.00706442}, {-4, 0.0120003}, {-3, -0.0367309}, {-2, -0.0488618}, {-1, 0.280931}, {0, 0.578163}, {1, 0.280931}, {2, -0.0488618}, {3, -0.0367309}, {4, 0.0120003}, {5, 0.00706442}, {6, -0.00274588}, {7, -0.00155701}, {8,0.000652922}, {9, 0.000361781}, {10, -0.000158601}, {11, -0.0000867523}, {12, 0.0000391428}, {13, 0.0000212211}, {14, -9.77137*10^-6}}

My code until now:

f[x_] := Piecewise[{{1, 0 <= x <= 1}}]

g[x_] := f[2 x] - f[2 x - 1];

MyWavelet[lim_Integer]["WaveletQ"] := True
MyWavelet[lim_]["OrthogonalQ"] := True

MyWavelet[lim_Integer]["PrimalLowpass", prec_ : MachinePrecision] := g[lim]

The Haar Wavelet in defined by it´s scaling function f(x) and has the wavelet function g(x). But obviously mathematica is waiting for some values that define the wavelet but how?

Answer 1

You almost have it working, just a couple minor fixes and it seems to work. Fist, in order to use DiscreteWaveletTransform I had to set it to be an orthogonal wavelet. Then I corrected your definition of the PrimalLowpass to include a second optional argument of the precision.

f[x_] := Piecewise[{{1, 0 <= x <= 1}}]
g[x_] := f[2 x] - f[2 x - 1]

MyWavelet[]["WaveletQ"] := True
MyWavelet[]["OrthogonalQ"] := True
MyWavelet[]["BiorthogonalQ"] := False
MyWavelet[]["WaveletFunction"] := g[#1] &
MyWavelet[]["FourierFactor"] := 1/Sqrt[2 Pi]
MyWavelet[]["FourierTransform"] := Sinc[#1/2] &
MyWavelet[]["PrimalLowpass", 
  prec_: MachinePrecision] := {{0, 0.5}, {1, 0.5}}

Now you can use DiscreteWaveletTransform,

WaveletScalogram[
 DiscreteWaveletTransform[Table[Sin[x^2], {x, -6, 6, 0.01}], 
  MyWavelet[]]]

As suggested by Sektor, you can test the validity of the transform by comparing the back-transformed data to the original data. This transform is good for lists of 63 data points or less when using ContinuousWaveletTransform,

ListPlot[
 Table[
  data = RandomReal[5, n];
  {n, Max[
    data - InverseContinuousWaveletTransform@
      ContinuousWaveletTransform[data, MyWavelet[]]]}, {n, 10, 100}]
 ]

Why it fails after that? I haven't the foggiest, but there is no such problem with the discrete transform.

ListPlot[
 Table[
  data = RandomReal[5, n];
  {n, Max[
    data - InverseWaveletTransform@
      DiscreteWaveletTransform[data, MyWavelet[]]]}, {n, 10, 1000}]
 ]

Answer 2

Until now I develop a partly solution of my own problem. With this solution someone can compute continuous wavelet transforms and inverse wavelet transforms as well as other examples below. But until now I dont fixed the problem with the low pass filter. As you can see I added the option for the filter coefficients but this is not working out. I think if this problem gets fixed someone can use every option of the framework with its own wavelet.

Some help to at least fix that would be realy nice.

f[x_] := Piecewise[{{1, 0 <= x <= 1}}]
g[x_] := f[2 x] - f[2 x - 1]

MyWavelet[]["WaveletQ"] := True
MyWavelet[]["OrthogonalQ"] := False
MyWavelet[]["BiorthogonalQ"] := False
MyWavelet[]["WaveletFunction"] := g[#1] &
MyWavelet[]["FourierFactor"] := 1/Sqrt[2 Pi]
MyWavelet[]["FourierTransform"] := Sinc[#1/2] &
MyWavelet[]["PrimalLowpass"] := {{0, 0.5}, {1, 0.5}}

With this options set someone can perform the following commands:

Plot[WaveletPsi[MyWavelet[], x], {x, -1, 2}, Exclusions -> None]

data = Table[Random[], {i, 7}]
CWTMyWavelet = ContinuousWaveletTransform[data, MyWavelet[]]
{0.607242, 0.866118, 0.898763, 0.802639, 0.445842,0.671422, 0.913055}

 InverseContinuousWaveletTransform[CWTMyWavelet]
{0.607242, 0.866118, 0.898763, 0.802639, 0.445842,0.671422, 0.913055}

WaveletScalogram[CWTMyWavelet, Ticks -> Full]

DiscreteWaveletTransform[data, MyWavelet]
DiscreteWaveletTransform::bdwave: The wavelet specification MyWavelet cannot be used to perform DiscreteWaveletTransform. >>