Merge pull request #76 from JaredW40/master

Updated Morse wavelet function and test functions

Author: BoundaryValueProblems

.github/workflows/CI.yml

@@ -10,7 +10,7 @@ on:
   pull_request: # run when pull requests are made
     paths:
       - src/*
-      - test/runtests.jl
+      - test/*
       - Manifest.toml
       - Project.toml
 jobs:
@@ -22,7 +22,7 @@ jobs:
       matrix:
         version:
           - "1.11"
-          - "1.10"
+          - "1.12"
         os:
           - ubuntu-latest
           - macOS-latest

docs/allInOne.jl

@@ -123,7 +123,7 @@ gr();
 Plots.reset_defaults()
 global_logger(Logging.SimpleLogger(stderr, Logging.Error))
 n = 2047;
-function mapTo(waveType, isReal = true, window = 1:2047; d = 1, γ = 4.0, β = 2.0, cf = 1.0, kwargs...)
+function mapTo(waveType, isReal = true, window = 1:2047; d = 1, γ = 2.0, β = 4.0, cf = 1.0, kwargs...)
     if waveType == Morse
         morse_wav = Morse(float(γ), float(β), float(cf))
         c = wavelet(morse_wav; kwargs...)
@@ -150,18 +150,18 @@ p2 = plot([real.(tmp) imag.(tmp)],
     	labels = ["real" "imaginary"],
     	ticks = nothing,
     	linewidth = 5)
-tmpMorse1 = mapTo(Morse, false; β=3, γ=10.0, cf=1.0, averagingLength=-1)[:, 2]
+tmpMorse1 = mapTo(Morse, false; γ=2.0, β=11.0, cf=1.0, averagingLength=-1)[:, 2]
 p3 = plot([real.(tmpMorse1) imag.(tmpMorse1)],
-    	title = "Morse (β=3, γ=10)",
+    	title = "Morse (γ=2, β=11)",
     	labels = ["real" "imaginary"],
     	ticks = nothing,
-    	linewidth = 4)
-tmpMorse2 = mapTo(Morse, false; β=1, γ=3.0, cf=1.0, averagingLength=-2)[:, 2] 
+    	linewidth = 5)
+tmpMorse2 = mapTo(Morse, false; γ=1, β=2, cf=1.0, averagingLength=-2)[:, 2] 
 p4 = plot([real.(tmpMorse2) imag.(tmpMorse2)],
-    	title = "Morse (β=1, γ=3)",
+    	title = "Morse (γ=1, β=2)",
     	labels = ["real" "imaginary"],
     	ticks = nothing,
-    	linewidth = 4)
+    	linewidth = 5)
 p5 = plot(mapTo(dog2; averagingLength = -1.5)[:, 2],
     	title = "derivative of gaussians (dog2)",
     	legend = false,

docs/mothers.svg

@@ -10,7 +10,7 @@ using ContinuousWavelets, Plots, Wavelets, FFTW, Logging
 using Plots; gr()
 Plots.reset_defaults()
 n= 2047;
-function mapTo(waveType, isReal = true, window = 1:2047; d = 1, γ = 4.0, β = 2.0, cf = 1.0, kwargs...)
+function mapTo(waveType, isReal = true, window = 1:2047; d = 1, γ = 2.0, β = 4.0, cf = 1.0, kwargs...)
     if waveType == Morse
         morse_wav = Morse(float(γ), float(β), float(cf))
         c = wavelet(morse_wav; kwargs...)
@@ -29,10 +29,10 @@ tmp = mapTo(Morlet(π), false; averagingLength = -0.2)[:, 2]
 p1 = plot([real.(tmp) imag.(tmp)], title = "Morlet", labels = ["real" "imaginary"], ticks = nothing, linewidth = 5)
 tmp = mapTo(paul2, false, averagingLength = -0.5)[:, 2]
 p2 = plot([real.(tmp) imag.(tmp)], title = "Paul 2", labels = ["real" "imaginary"], ticks = nothing, linewidth = 5)
-tmpMorse1 = mapTo(Morse, false; β=3, γ=10.0, cf=1.0, averagingLength=-1)[:, 2]
-p3 = plot([real.(tmpMorse1) imag.(tmpMorse1)], title = "Morse (β=3, γ=10)", labels = ["real" "imaginary"], ticks = nothing, linewidth = 4)
-tmpMorse2 = mapTo(Morse, false; β=1, γ=3.0, cf=1.0, averagingLength=-2)[:, 2] 
-p4 = plot([real.(tmpMorse2) imag.(tmpMorse2)], title = "Morse (β=1, γ=3)", labels = ["real" "imaginary"], ticks = nothing, linewidth = 4)
+tmpMorse1 = mapTo(Morse, false; γ=2.0, β=11.0, cf=1.0, averagingLength=-1)[:, 2]
+p3 = plot([real.(tmpMorse1) imag.(tmpMorse1)], title = "Morse (γ=2, β=11)", labels = ["real" "imaginary"], ticks = nothing, linewidth = 4)
+tmpMorse2 = mapTo(Morse, false; γ=1.0, β=2.0, cf=1.0, averagingLength=-2)[:, 2] 
+p4 = plot([real.(tmpMorse2) imag.(tmpMorse2)], title = "Morse (γ=1, β=2)", labels = ["real" "imaginary"], ticks = nothing, linewidth = 4)
 p5 = plot(mapTo(dog2; averagingLength = -1.5)[:, 2], title = "derivative of gaussians (dog2)", legend = false, ticks = nothing, linewidth = 5)
 p6 = plot(mapTo(cHaar, true; averagingLength = 1)[:, 2], title = "Haar", legend = false, ticks = nothing, linewidth = 5)
 p7 = plot(mapTo(cBeyl, true; d = 1, averagingLength = -0)[:, 2], title = "Beylkyin", legend = false, ticks = nothing, linewidth = 5)
@@ -46,5 +46,5 @@ plot(p1, p2, p3, p4, p5, p6, p7, p8, p9, p10, p11, p12, layout = (3, 4), size =
 savefig("mothers.svg")#hide
 ```
 ![](mothers.svg)
-Above are examples of every mother wavelet family defined in this package; the only analytic and/or complex wavelets are the `Morlet` and the `Paul` wavelet families.
+Above are examples of every mother wavelet family defined in this package; the only analytic and/or complex wavelets are the `Morlet`, `Paul` and `Morse` wavelet families.
 Once you have chosen a type of wavelet, this is used to construct the more specific CWT, which specifies more details of the transform, such as frequency spacing, whether to include an averaging filter or not, a frame upper bound, etc.

docs/src/coreType.md

@@ -69,7 +69,7 @@ are:
     compactness parameter, and `cf` is a parameter that determines the 
     frequency at which the magnitude of the wavelet is maximized. 
 
-    ``\\psi\\hat(\\omega) \\propto \\chi_{\\omega \\geq 0}(\\omega) \\omega^\\textrm{be}\\textrm{e}^{-(2\\pi\\omega)^{\\textrm{ga}}}\\textrm{, } be, ga \\geq 0``
+    ``\\psi\\hat(\\omega) \\propto \\chi_{\\omega \\geq 0}(\\omega) \\omega^\\textrm{be}\\textrm{e}^{-\\omega^{\\textrm{ga}}}\\textrm{, } ga, be \\geq 0``
 
 - `Dog{N}`: Derivative of a Gaussian, where N is the number of
     derivatives. `dogn` for `n` in `0:6`. The Sombrero/mexican hat/Marr wavelet

src/ContinuousWavelets.jl

@@ -48,27 +48,29 @@ function mother(this::CWT{W,T,Morse,N},
 
     ga = this.waveType.ga
     be = this.waveType.be
-    cf = this.waveType.cf
+    # cf = this.waveType.cf
     p = this.p
 
     fo = morsefreq(this)
-    fact = cf / fo
+    # fact = cf / fo
 
     #  ω = LinRange(0,1-(1/len),len)
     # om = 2 * pi * ω./ fact / max(1, s)
-    #om = 2 * pi * (ω / s)./ fact
+    # om = 2 * pi * (ω / s)./ fact
     # om = (ω / s) / cf
     # om = (ω / s) / fact
 
     om = ω / s
+    om_safe = max.(om, eps())
 
     if be == 0
         daughter = @. 2 * exp(-om^ga)
     else
-        daughter = @. 2 * exp(-be * log(fo) + fo^ga + be * log(om) - om^ga)
+        daughter = @. 2 * (om_safe / fo)^be * exp(-(om^ga - fo^ga))
     end
-
-    daughter[1] = 1 / 2 * daughter[1] # Due to unit step function
+    daughter[1] = 0
+    
+    # daughter[1] = 1 / 2 * daughter[1] # Due to unit step function
     # Ensure nice lowpass filters for beta=0;
     # Otherwise, doesn't matter since wavelets vanishes at zero frequency
 

src/createWavelets.jl

@@ -81,8 +81,7 @@ getMinScaling(c::CWT{W,T,M,N}) where {W,T,N,M} = 0 # by default all scales are a
 getMinScaling(c::CWT{W,T,<:Morlet,N}) where {W,T,N} = 1 / (c.β)^0.8 # morlet is slightly too large at the boundaries by default
 getMinScaling(c::CWT{W,T,<:Paul,N}) where {W,T,N} = 2 / (2c.α + 1) # Paul presents some difficulties, as the decay changes quickly (like 1/t^(α+1))
 getMinScaling(c::CWT{W,T,<:Dog,N}) where {W,T,N} = 2 # like morlet, the decay for Dog is exponential and consistent across derivatives
-# getMinScaling(c::CWT{W,T,<:Morse,N}) where {W,T,N,M} = log2(2 * morsefreq(c))
-getMinScaling(c::CWT{W,T,<:Morse,N}) where {W,T,N} = 5 * morsefreq(c)
+getMinScaling(c::CWT{W,T,<:Morse,N}) where {W,T,N} = log2(2 * morsefreq(c))
 
 function varianceAdjust(this::CWT{W,T,M,N}, totalWavelets, nOct) where {W,T,N,M}
     # increases the width of the wavelets by σ[i] = (1+a(total-i)ᴾ)σₛ

src/utils.jl

@@ -37,7 +37,7 @@ function Morse_convert(ga::Real, be::Real, cf::Real)
     Morse(ga, be, cf)
 end
 
-Morse() = Morse_convert(4, 2, 1)
+Morse() = Morse_convert(2, 4, 1)
 class(::Morse) = "Morse";
 name(::Morse) = "morse";
 vanishingmoments(::Morse) = 0;

src/waveletTypes.jl

@@ -38,19 +38,25 @@
                 end-floor(Int, n / 2):end],
             :]
         # make sure the wavelets don't have significant support out into the mirrored signal
-        if !(bc isa PerBoundary) && !(wave isa Paul) && size(Ŵ, 2) > 1
+        if !(bc isa PerBoundary) && !(wave isa Paul) && !(wave isa Morse) && size(Ŵ, 2) > 1
             @test max([norm(nonSupported[:, i], Inf) / norm(supported[:, i], Inf) for
                        i = 2:size(spaceWaves, 2)]...) ≤ 1e-2
-        elseif wave isa Paul && size(Ŵ, 2) > 1 # the Paul wavelet decay is just too slow to avoid bleedover, so best to ignore the averaging size
+        elseif !(bc isa PerBoundary) && (wave isa Paul || wave isa Morse) && 
+                size(Ŵ, 2) > 1 && ave > 0
             @test max([norm(nonSupported[:, i], Inf) / norm(supported[:, i], Inf) for
                        i = 2:size(spaceWaves, 2)]...) ≤ 1e-1
         end
         # if the averaging length is 0 and its only averaging at the sizes given, something is wrong
         @test size(Ŵ, 2) > 1 || ave > 0
+
         # make sure that the highest frequency support is small relative to the space domain
-        # Guaranteed in a very different way for ContOrtho, no guarantees made for just averaging, and the averaging length needs to leave 6.5 octaves
-        if !(wave isa ContOrtho) && size(Ŵ, 2) > 1 && log2(n) - ave > 6.5
-            @test max(abs.(Ŵ[end, :])...) / norm(supported, Inf) ≤ 1e-1
+        # Guaranteed in a very different way for ContOrtho, no guarantees made for just averaging,
+        # and the averaging length needs to leave 6.5 octaves.
+        # The Morse wavelet with γ < 3 has sub-Gaussian spectral decay and cannot be guaranteed
+        # to vanish at Nyquist, similar to the Paul wavelet.
+        if !(wave isa ContOrtho) && size(Ŵ, 2) > 1 && log2(n) - ave > 6.5 &&
+                !(wave isa Paul) && !(wave isa Morse && wave.ga < 3)
+            @test max(abs.(Ŵ[end, :])...) / norm(supported, Inf) ≤ 1e-1
         end
 
         # check that the father wavelet is actually positive