Skip to main content

Continuous wavelet transform with complex Morlet function


The complex Morlet function is defined as:


$$Ψ(t,f_c,f_b)= \frac{1}{\sqrt[]{ \pi f_{b} } }\exp(-t^2/f_b)\exp(\jmath 2πf_ct)$$


where $f_b$ and $f_c$ are two important parameters in modifying the complex Morlet wavelet. It seems that Mathematica doesn't support complex Morlet transform and Its only support real morlet function that I am not interested to use. I'm into complex wavelet function. Mathematica only has Gabor transform for complex wavelets, and Gabor transform only has one parameter to be tuned.
so I need complex morlet function to run continues wavelet transform. Also I want to define $f_b$ and $f_c$ of the complex morlet function myself.
Can I make a complex Morlet wavalet transform by changing Gabor's parameter? How can I change $f_b$ and $f_c$ in it?
can I define a new wavelet exactly like the equation of complex morlet?


P.S: Actually I am a MATLAB user and as such I don't really know anything about the flexibility of Mathematica, but the reason why I came here is that Mathematica has the InverseContinuousWaveletTransform.



Answer




EDIT:


First, a note: As the usage of options, parameters and functions listed below is not documented, be advised that they still need proper tuning and/or may not work at all.


CMorletWavelet[]["WaveletQ"] := True
CMorletWavelet[]["OrthogonalQ"] := False
CMorletWavelet[]["BiorthogonalQ"] := False
CMorletWavelet[]["WaveletFunction"] := 1/Sqrt[π] Exp[2 I π 2 #1] Exp[-#1^2] &
CMorletWavelet[]["FourierFactor"] := 4 π/(6 + Sqrt[2 + 6^2])
CMorletWavelet[]["FourierTransform"] := Function[{Wavelets`NonOrthogonalWaveletsDump`wt, 
   Wavelets`NonOrthogonalWaveletsDump`s},
  π^(-1/4)HeavisideTheta[Wavelets`NonOrthogonalWaveletsDump`wt + $MachineEpsilon] 

  Exp[-(1/2) (Wavelets`NonOrthogonalWaveletsDump`wt Wavelets`NonOrthogonalWaveletsDump`s
  - π Sqrt[2/Log[2]])^2]]

Now you can use the built-in wavelet-related functions:


Plot[{Re@WaveletPsi[CMorletWavelet[], x], Im@WaveletPsi[CMorletWavelet[], x]},
     {x, -5, 5}, PlotRange -> All, Frame -> True, GridLines -> Automatic, 
     PlotStyle -> {Blue, {Red, Dashed}}]


wavelet_plot




snd = Play[Sum[Sin[2000 2^t n t], {n,5 }], {t, 2, 3}]

csd = ContinuousWaveletTransform[snd, CMorletWavelet[]]

WaveletScalogram[csd]


Sound scalogram




InverseContinuousWaveletTransform[csd, CMorletWavelet[]]


Play sound



This sound compression works just fine !


(* A simple example *)
cwd = ContinuousWaveletTransform[Range[10], CMorletWavelet[]]
WaveletScalogram[cwd]



Scalogram



 InverseContinuousWaveletTransform[cwd, CMorletWavelet[]]


{1., 2., 3., 4., 5., 6., 7., 8., 9., 10.}

This works just as expected, but using numbers larger than 63 results in ..


 cwd = ContinuousWaveletTransform[Range[64], CMorletWavelet[]]

 WaveletScalogram[cwd]


Scalogram



 InverseContinuousWaveletTransform[cwd, CMorletWavelet[]]


{0.500005, 4.38214, 6.69958, 10.625, 12.6907, 16.5033, 18.2989, 
21.8762, 23.3564, 26.6196, 27.7395, 30.6377, 31.3658, 33.8706, 

34.1929, 36.2965, 36.2168, 37.9296, 37.4675, 38.8152, 38.0038, 
39.0243, 37.9069, 38.647, 37.274, 37.7859, 36.2116, 36.551, 34.8323, 
35.0564, 33.2508, 33.4173, 31.5827, 31.7492, 29.9436, 30.1677, 
28.449, 28.7884, 27.2141, 27.726, 26.353, 27.0931, 25.9757, 26.9962, 
26.1848, 27.5325, 27.0704, 28.7832, 28.7035, 30.8071, 31.1294, 
33.6342, 34.3623, 37.2605, 38.3804, 41.6436, 43.1238, 46.7011, 
48.4967, 52.3093, 54.375, 58.3004, 60.6179, 64.5}

One of the reasons for this lies in the fact that I used the Fourier Transform of the original MorletWavelet which is a built-in predicate and quite different in implementation from the one I used. There are probably other parameters I need to set up properly, but I can't seem to find them, because, like I said, the usage is undocumented.




I know you came here because of the InverseContinuousWaveletTransform, but at that time of the day, or should I say night, I can't really think any more and will continue when I have more time to do so, unfortunately...


Note: As you are a MATLAB user I implemented the Complex Morlet wavelet according to THEIR documentation.


Preliminaries


For simplicity we assume that smallest wavelet scale is equal to 1 and we use a rather short data set.


I also used the following pages from the documentation (A-Z)



Implementation


(* Example data set *)

data = {1, 2, 3, 4};


(* Parameters *)

noct = Floor@Log[2, (data // Length)/2]


1

nvoc = 4;


(* Scaling parameter *)

s[oct_, voc_] := N[2^(oct - 1) 2^(voc/nvoc)]

(* Defining the wavelet function *)

ComplexMorlet[n_, band_, centerFreq_] := 
    1/Sqrt[π band] Exp[2 I π centerFreq n] Exp[-n^2/band]

(* Example expansion *)


ComplexMorlet[x, 1, 2]


E^(4 I π x - x^2)/Sqrt[π]

Plot[{Re@ComplexMorlet[x, 1, 2], Im@ComplexMorlet[x, 1, 2]}, {x, -3, 3},
     PlotStyle -> {Blue, {Red, Dashed}}, PlotRange -> All, 
     Frame -> True, GridLines -> Automatic]



Plot



(* Wavelet transform of a sampled sequence *)

 w[u_, oct_, voc_] := 1/s[oct, voc] Sum[data[[k]]
     Conjugate[ComplexMorlet[(k - u)/s[oct, voc], 1, 2]], {k, 1, data // Length}]

(* Performing the wavelet transform on our example data set *)


Table[w[k, 1, voc], {k, data // Length}, {voc, 4}]


{{0.228074 + 0.361025 I, 0.0610598 - 0.123408 I, 
     0.283659 - 0.583475 I, 1.15175 + 3.47516*10^-16 I},
   {0.486587 + 0.340747 I, 0.0693978 - 0.058132 I, 0.786587 - 0.662852 I, 
     1.85808 + 3.10964*10^-16 I}, 
   {0.821662 + 0.446737 I, -0.0236108 - 0.295969 I, 1.47435 - 0.380752 I, 
     2.26824 + 5.67838*10^-17 I}, 
   {1.57014 - 0.595682 I, 1.02407 + 0.281895 I, 1.47482 + 0.762858 I, 

     2.02475 - 2.84949*10^-16 I}}

(* Wavelet Scalogram using ComplexMorlet[x, 1, 2] *)

WaveletScalogram@ContinuousWaveletData[
 {{1, 1} -> {0.22807383843702972` + 0.36102529036876024` I, 
       0.06105984372279422` - 0.12340783119864777` I, 
       0.28365883675526904` - 0.5834746966816698` I, 
       1.1517469935306757` + 3.4751640646106677`*^-16 I},
  {1, 2} -> {0.4865866432814967` + 0.3407467247569226` I, 

       0.06939782717412021` - 0.05813200432524761` I, 
       0.7865874222126943` - 0.6628516103818837` I, 
       1.8580796599037956` + 3.1096385445125467`*^-16 I},
  {1, 3} -> {0.8216617511105463` + 
       0.44673675942817265` I, -0.02361080340458542` - 
       0.2959689122870983` I, 
       1.4743517412825382` - 0.3807516306374966` I, 
       2.26823511807995` + 5.678382044215492`*^-17 I},
  {1, 4} -> {1.570143054029254` - 0.5956822545417808` I, 
       1.024067417876664` + 0.2818946441776095` I, 

       1.4748223337693926` + 0.7628582023394818` I, 
       2.024752422313301` - 2.849488941725102`*^-16 I}}]


WaveletScalogram



(* Wavelet Scalogram using ComplexMorlet[x, 1, 10] *)

WaveletScalogram@ContinuousWaveletData@
 {{1, 1} -> {0.11634486079523618` - 0.17990847470866217` I, 

       0.9410569485064904` - 0.3524175549056541` I, 
       0.9995892268140318` + 0.3575695443712028` I, 
       1.1517469935306757` + 2.5826325630023094`*^-15 I},
 {1, 2} -> {0.2085276338912312` - 0.15114828701865127` I, 
      1.8062819251440743` - 0.3772206439472593` I, 
      1.813592761954768` + 0.36136020250254647` I, 
      1.8580796599037956` + 1.5548192722562736`*^-15 I},
 {1, 3} -> {0.2547509048762912` - 0.27877696228455096` I, 
      2.5401537117071564` - 0.16692666476822` I, 
      2.402824979378204` + 0.10553538050034861` I, 

      2.26823511807995` + 2.8391910221077465`*^-16 I},
 {1, 4} -> {1.3309683457126755` + 0.3296339838999044` I, 
      2.319228847343012` + 0.4019097092762081` I, 
      2.1426745757435186` - 0.3492240227193354` I, 
      2.024752422313301` - 1.6360071035367952`*^-15 I}}


WaveletScalogram 2



Comments

Post a Comment

Popular posts from this blog

plotting - Plot 4D data with color as 4th dimension

I have a list of 4D data (x position, y position, amplitude, wavelength). I want to plot x, y, and amplitude on a 3D plot and have the color of the points correspond to the wavelength. I have seen many examples using functions to define color but my wavelength cannot be expressed by an analytic function. Is there a simple way to do this? Answer Here a another possible way to visualize 4D data: data = Flatten[Table[{x, y, x^2 + y^2, Sin[x - y]}, {x, -Pi, Pi,Pi/10}, {y,-Pi,Pi, Pi/10}], 1]; You can use the function Point along with VertexColors . Now the points are places using the first three elements and the color is determined by the fourth. In this case I used Hue, but you can use whatever you prefer. Graphics3D[ Point[data[[All, 1 ;; 3]], VertexColors -> Hue /@ data[[All, 4]]], Axes -> True, BoxRatios -> {1, 1, 1/GoldenRatio}]
Post a Comment
Read more

plotting - Filling between two spheres in SphericalPlot3D

Manipulate[ SphericalPlot3D[{1, 2 - n}, {θ, 0, Pi}, {ϕ, 0, 1.5 Pi}, Mesh -> None, PlotPoints -> 15, PlotRange -> {-2.2, 2.2}], {n, 0, 1}] I cant' seem to be able to make a filling between two spheres. I've already tried the obvious Filling -> {1 -> {2}} but Mathematica doesn't seem to like that option. Is there any easy way around this or ... Answer There is no built-in filling in SphericalPlot3D . One option is to use ParametricPlot3D to draw the surfaces between the two shells: Manipulate[ Show[SphericalPlot3D[{1, 2 - n}, {θ, 0, Pi}, {ϕ, 0, 1.5 Pi}, PlotPoints -> 15, PlotRange -> {-2.2, 2.2}], ParametricPlot3D[{ r {Sin[t] Cos[1.5 Pi], Sin[t] Sin[1.5 Pi], Cos[t]}, r {Sin[t] Cos[0 Pi], Sin[t] Sin[0 Pi], Cos[t]}}, {r, 1, 2 - n}, {t, 0, Pi}, PlotStyle -> Yellow, Mesh -> {2, 15}]], {n, 0, 1}]
Post a Comment
Read more

plotting - Mathematica: 3D plot based on combined 2D graphs

I have several sigmoidal fits to 3 different datasets, with mean fit predictions plus the 95% confidence limits (not symmetrical around the mean) and the actual data. I would now like to show these different 2D plots projected in 3D as in but then using proper perspective. In the link here they give some solutions to combine the plots using isometric perspective, but I would like to use proper 3 point perspective. Any thoughts? Also any way to show the mean points per time point for each series plus or minus the standard error on the mean would be cool too, either using points+vertical bars, or using spheres plus tubes. Below are some test data and the fit function I am using. Note that I am working on a logit(proportion) scale and that the final vertical scale is Log10(percentage). (* some test data *) data = Table[Null, {i, 4}]; data[[1]] = {{1, -5.8}, {2, -5.4}, {3, -0.8}, {4, -0.2}, {5, 4.6}, {1, -6.4}, {2, -5.6}, {3, -0.7}, {4, 0.04}, {5, 1.0}, {1, -6.8}, {2, -4.7}, {3, -1....
Post a Comment
Read more