Cone of influence

[henry2004y]

Hello,

I am very excited to see the scaleogram demos from continuous wavelet transform!

One step further, I would love to see the addition of the cone of influence on top of the scaleogram. When searching through the source codes, I found one function caveats which returns coi as the 3rd argument:

ContinuousWavelets.jl/src/apply.jl

Line 233 in 6e4b54d

function caveats(n1, c::CWT{W}; J1::Int64=-1, dt::S=1/1000, s0::V=NaN) where {S<:Real, W<:WaveletBoundary, V <: Real}

However, there is currently no test or example of calling this function.

[dsweber2]

It was part of a previous implementation which output that along with the actual cwt. I'll bring it up to date, as it is actually useful.

[dsweber2]

This is now updated in bcaca5b so that the cone of influence is included. It's a pretty slow implementation that directly computes it from the given wavelets, rather than using the explicit constructions in, e.g. Torrence and Compo, but it will give exact boundaries.

[henry2004y]

Thanks for separating this piece from the whole caveat calculation!

I've been testing on this cone of influence. The current return value is of type BitArray, which requires some extra work to actually plot the COI on top of the spectrum. Here I use PyPlot.jl as an example.

using PyPlot, Wavelets, ContinuousWavelets

n = 2047
t = range(0, n/1000, length=n) # 1kHz sampling rate
f = testfunction(n, "Doppler")

fig, ax = subplots(3,1, constrained_layout=true)

p1 = ax[1].plot(t, f)
ax[1].set_title("Doppler")

c = wavelet(Morlet(π), β=2)
res = ContinuousWavelets.cwt(f, c)

freqs = getMeanFreq(ContinuousWavelets.computeWavelets(n, c)[1])
freqs[1] = 0

p1 = ax[2].pcolormesh(t, freqs, abs.(res)', cmap = matplotlib.cm.turbo, shading="nearest")

# 2047*33
coiMatrix = ContinuousWavelets.directCoiComputation(n, c; coiTolerance=exp(-1))

function extract_coi_line(coiMatrix, freqs)
   coi_index = ones(Int64, size(coiMatrix, 1))

   for i in axes(coiMatrix,1)
      for j in Iterators.reverse(axes(coiMatrix,2))
         if coiMatrix[i,j] == 1
            coi_index[i] = j
            break
         end
      end
   end

   freqs[coi_index]
end

coi = extract_coi_line(coiMatrix, freqs)

ax[2].fill_between(t, coi, y2=0; color="white", alpha=0.7)

p2 = ax[3].pcolormesh(t, freqs, coiMatrix', shading="nearest")

ax[2].set_xlabel("time (s)")
ax[2].set_ylabel("frequency (Hz)")

colorbar(p1; ax=ax[2])
colorbar(p2; ax=ax[3])

which generates

I hope to add something similar to this extract_coi_line to the package. Do you think it's proper?

---

Another tiny issue I've noticed from the above demo plot is that for the lowest frequency coefficients, coiMatrix returns 1s but not 0s, which creates a stripe in the middle at the bottom. Do you think there's an error in the calculation?

---

I've also noticed you had something not complete in testing the COI. Maybe we can simply add a test for this based on the envelop of the cone of influence?

[dsweber2]

I hope to add something similar to this extract_coi_line to the package. Do you think it's proper?

Yes, that would probably be helpful, thanks (though camel case rather than snake case, and probably output the indices too). I'm thinking about adding an actual scalogram function using Plots recipes, and this would be useful for that.

Another tiny issue I've noticed from the above demo plot is that for the lowest frequency coefficients, coiMatrix returns 1s but not 0s, which creates a stripe in the middle at the bottom. Do you think there's an error in the calculation?

Not sure I quite see what you're getting at. If you're concerned about the fact that the width decreases between the first and second wavelet, the first wavelet is the father wavelet, which may have narrower time support than the lowest frequency wavelets, as is happening in this case.

I've also noticed you had something not complete in testing the COI. Maybe we can simply add a test for this based on the envelop of the cone of influence?

Oh, I hadn't intended to actually commit those and should remove them. I was planning on doing a comparison between the theoretical approach and this direct approach, though other tests would be good too.

[henry2004y]

Not sure I quite see what you're getting at. If you're concerned about the fact that the width decreases between the first and second wavelet, the first wavelet is the father wavelet, which may have narrower time support than the lowest frequency wavelets, as is happening in this case.

No, that's not what I meant. If you look at the 3rd subplot, that's the raw BitArray display of the cone of influence, with 0s meaning "this probably can be trusted" and 1s "no this is doubtful". Between about 0.1s - 0.3s and 1.7s - 2.0s at the lowest y row, which is also the lowest frequency, they are 0s, while at the 2nd lowest y row they are 1s, which is very confusing to me.
Maybe that's related to the warning earlier?

julia> res = ContinuousWavelets.cwt(f, c)
┌ Warning: the lowest frequency wavelet has more than 1% its max at zero, so it may not be analytic. Think carefully
│   lowAprxAnalyt = 0.06186323501016359
└ @ ContinuousWavelets ~/.julia/packages/ContinuousWavelets/lDZFX/src/sanityChecks.jl:6

Or is it just the DC component, which shouldn't be plotted at all?

I'm thinking about adding an actual scalogram function using Plots recipes, and this would be useful for that.

One issue of adding a plot recipe is that we have all the pieces, but they are bundled together. The time range t, frequency range freqs, and CWT coefficients res, are all regular data types. For a recipe, we need something like

struct spectra
  time
  freq
  coef
end

and then

@recipe function f(myspectra::spectra, ...; ...) end

Also, we don't want to do computations inside the recipes, so the WT shall be done beforehand.

[dsweber2]

Or is it just the DC component, which shouldn't be plotted at all?

This. The DC component/father/averaging wavelet for this particular setup also has pretty wide frequency support (so its narrow in space)