Merge pull request #54 from dsweber2/fixingMeanFreq

mean power (rather than 1-norm) and linear interpolation
UCD4IDS · Jun 3, 2024 · 4b92937 · 4b92937 · dsweber2 · Jun 3, 2024
2 parents dbb242b + 92a1269
commit 4b92937
Show file tree

Hide file tree

Showing 4 changed files with 11 additions and 10 deletions.
diff --git a/Project.toml b/Project.toml
@@ -1,10 +1,11 @@
 name = "ContinuousWavelets"
 uuid = "96eb917e-2868-4417-9cb6-27e7ff17528f"
 authors = ["dsweber2 <david.weber2@protonmail.com> and contributors"]
-version = "1.1.4"
+version = "1.1.5"
 
 [deps]
 AbstractFFTs = "621f4979-c628-5d54-868e-fcf4e3e8185c"
+DSP = "717857b8-e6f2-59f4-9121-6e50c889abd2"
 Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4"
 FFTW = "7a1cc6ca-52ef-59f5-83cd-3a7055c09341"
 Interpolations = "a98d9a8b-a2ab-59e6-89dd-64a1c18fca59"

diff --git a/src/createWavelets.jl b/src/createWavelets.jl
@@ -186,7 +186,7 @@ function computeWavelets(n1::Integer,
     isAve = !(typeof(c.averagingType) <: NoAve)
     # I guess matlab did occasionally do something useful
 
-    ω = computeOmega(n1, nSpace, n)
+    ω = computeOmega(n1, n)
     daughters = analyticOrNot(c, n, totalWavelets)
 
     # if the nOctaves is small enough there are none not covered by the
@@ -230,7 +230,7 @@ function computeWavelets(n1::Integer,
     # indicates whether we should keep a spot for the father wavelet
     isAve = !(typeof(c.averagingType) <: NoAve)
 
-    ω = computeOmega(n1, nSpace, n)
+    ω = computeOmega(n1, n)
     daughters = zeros(nSpace, totalWavelets)
     φ, ψ, ψLen = getContWaveFromOrtho(c, nSpace)
     itpψ = genInterp(ψ)
@@ -287,7 +287,7 @@ function analyticOrNot(c::CWT{W,T,<:Union{Morlet,Paul,Morse},N},
     return daughters
 end
 
-function computeOmega(nOriginal, nSpace, nFreq)
+function computeOmega(nOriginal, nFreq)
     range(0, nOriginal >> 1 + 1, length = nFreq) # max size is the last frequency in the rfft of the original data size
 end
 # convenience methods

diff --git a/src/utils.jl b/src/utils.jl
@@ -247,13 +247,13 @@ end
 """
     getMeanFreq(Ŵ, fsample=2000) -> arrayOfFreqs
 
-Calculate each of the mean frequencies of a collection of analytic or real wavelets `Ŵ`.
+Calculate each of the mean frequencies of a collection of analytic or real wavelets `Ŵ`, using the power spectral density.
 The default sampling rate `fsample=2kHz`, so the maximum frequency is 1kHz.
 """
 function getMeanFreq(Ŵ::Array, fsample = 2000)
-    eachNorm = [norm(w, 1) for w in eachcol(Ŵ)]
+    eachNorm = [norm(w, 2)^2 for w in eachcol(Ŵ)]
     freqs = range(0, fsample / 2, length = size(Ŵ, 1))
-    return map(ŵ -> sum(abs.(ŵ) .* freqs), eachcol(Ŵ)) ./ eachNorm
+    return map(ŵ -> sum(abs2.(ŵ) .* freqs), eachcol(Ŵ)) ./ eachNorm
 end
 
 function getMeanFreq(n1, cw::CWT, fsample = 2000)
@@ -306,8 +306,8 @@ end
 
 
 # create interpolater for the orthogonal cases
-genInterp(ψ) = interpolate(ψ, BSpline(Quadratic(Reflect(OnGrid()))))
-
+#genInterp(ψ) = scale(interpolate(ψ, BSpline(Cubic(Reflect(OnGrid())))), 1:length(ψ))
+genInterp(ψ) = linear_interpolation(1:length(ψ),ψ)
 
 
 

diff --git a/src/waveletTypes.jl b/src/waveletTypes.jl
@@ -1,5 +1,5 @@
 struct Morlet <: ContWaveClass
-    σ::Float64 # σ is the time/space trade-off. as σ->0, the spacial resolution increases; below 5, there is a danger of being non-analytic. Default is 5.8
+    σ::Float64 # σ is the time/space trade-off. as σ->0, the spacial resolution increases; below 5, there is a danger of being non-analytic. Default is 2π
     κσ::Float64
     cσ::Float64
 end