4

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?

Dan
  • 233
  • 1
  • 8
  • 1
    "I tried to implement..." - so why not include the code you wrote, which might help us diagnose the problem? – J. M.'s missing motivation Mar 03 '16 at 11:14
  • To come up with an answer, most of us (and certainly me), would need to experiment with your code. But you don't give us any code. You are unlikely to get an answer unless you edit your question to include *all* the code needed to reproduce your problem. – m_goldberg Mar 03 '16 at 11:20
  • I added some code. I hope u can open the topic again. ty – Dan Mar 04 '16 at 12:59
  • It might be useful, http://mathematica.stackexchange.com/questions/33549/continuous-wavelet-transform-with-complex-morlet-function – s.s.o Mar 04 '16 at 14:34
  • 1
    @Dan, In every example on that page, the wavelet primallowpass returns a list, but yours just returns a number. – Jason B. Mar 04 '16 at 14:59
  • Looked just a bit farther into it, I don't understand how they are constructed from the g[x] function you have there, but if you look at this page you see that you need to set up the low-pass filter as {{0, 0.5}, {1, 0.5}}. The Haar wavelet as usually defined doesn't need an input parameter. – Jason B. Mar 04 '16 at 15:42
  • Also, on the tutorial you linked it states "Both properties PrimalLowPass and DualLowPass are expected to return a list of the form....." – Jason B. Mar 04 '16 at 15:43

3 Answers3

5

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[]]]

enter image description here

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}]
 ]

enter image description here

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}]
 ]

enter image description here

Jason B.
  • 68,381
  • 3
  • 139
  • 286
3

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]

enter image description here

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]

enter image description here

DiscreteWaveletTransform[data, MyWavelet]
DiscreteWaveletTransform::bdwave: The wavelet specification MyWavelet cannot be used to perform DiscreteWaveletTransform. >>
Dan
  • 233
  • 1
  • 8
0

In the definition of Haar wavelet transformation, the definition of Fourier factor and Fourier transform, I think it should be in the range of -1/2 <= x <= 1/2 not for 0 <= x <= 1. I mean, according to MyWavelet[]["FourierFactor"] := 1/Sqrt[2 Pi] and MyWavelet[]["FourierTransform"] := Sinc[#1/2], f[x_] should be like this:

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

not in the form of

f[x_] := Piecewise[{{1, 0 <= x <= 1}}]
creidhne
  • 5,055
  • 4
  • 20
  • 28