Feature request: cross wavelet tranform

[henry2004y]

Hi,

Thanks for sharing this package!

I am currently learning about wavelets and their practical usage in Julia. Based I what I've learned so far, the wavelet related functionalities currently in Julia are limited compared to more mature toolboxes like in MATLAB.

One thing I'm wondering specifically is the cross wavelet transform, or wavelet coherence described here. Is it possible to quickly implement a similar method within this package?

[dsweber2]

the wavelet related functionalities currently in Julia are limited compared to more mature toolboxes like in MATLAB.

the pros and cons of working in new open source languages.

Is it possible to quickly implement a similar method within this package?

It should be, as it's just a smoothed version of the product of the wavelets. For each wavelet mother there's an associated averaging father wavelet which would be a natural candidate for the smoothing. A pull request implementing that would be welcome if you want the practice, otherwise I'll get around to it sometime this week.

[henry2004y]

I will try to understand what's already there. Can I switch from Fourier analysis to wavelet analysis in a week :thinking:?

[dsweber2]

If you want to do it, I can certainly wait for that. The implementation shouldn't be too bad, as its just

The wavelet cross spectrum of two time series, x and y, is:

Cxy(a,b)=S(C^*x(a,b)Cy(a,b))

Cx(a,b) and Cy(a,b) denote the continuous wavelet transforms of x and y at scales a and positions b. The superscript * is the complex conjugate, and S is a smoothing operator in time and scale.

And for the Coherence

|Cxy(a,b)|^2 / S((C^*x(a,b)) / S(|Cy(a,b)|^2)

from the matlab documentation. It would be nice if github did latex.

[henry2004y]

I am not sure about what to use for the smoothing step. I would envision something like

using ContinuousWavelets, Wavelets

# time-series
f1 = testfunction(n, "Doppler")
f2 = copy(f1)

c = wavelet(Morlet(π), β=2)

c1 = ContinuousWavelets.cwt(f1, c)
c2 = ContinuousWavelets.cwt(f2, c)

# coherence
# smooth(c1'.*c2) .^2 ./ ( smooth(c1.^2) .* smooth(c2.^2) )
smooth(conj(c1).*c2) .^2 ./ ( smooth(c1.^2) .* smooth(c2.^2) )

should do the trick. I think it may be better for you @dsweber2 to handle the implementation 😅 .

[henry2004y]

And speaking of LaTeX support in github markdown, there are some tools like md-math.netliify.app for generating formulas for rendering in html.

It is still not optimal, but better than nothing :)

[henry2004y]

Now I'm trying to understand one step further. When we construct the continuous wavelet in the package, there is a property called averagingType which can be Father, Dirac, or NoAve. In createWavelets.jl, a father() method is provided for creating the averaging function, which I suppose can be used as the smoothing operator in calculating the wavelet coherence.

Here comes the part that is slightly confusing to me. When performing the transform in apply.jl, the averagingType is also passed as an argument to analyticTransformReal!, analyticTransformComplex!, and otherwiseTransform!,

https://github.com/UCD4IDS/ContinuousWavelets.jl/blob/9aceade7808f3b8748cc00c86679680f51192496/src/apply.jl#L131

but it is not actually used inside those functions. Instead, when the averagingType is either Father or Dirac, the same averaging process is performed

https://github.com/UCD4IDS/ContinuousWavelets.jl/blob/9aceade7808f3b8748cc00c86679680f51192496/src/apply.jl#L140

So my questions are:

1. Is averaging already done when doing CWT if averagingType is either Father or Dirac?
2. It seems like averagingType is simply used for multiple-dispatch in functions like analyticTrasnformReal!. In this case, maybe we can just keep the type in the arguments, something like

analyticTransformComplex!(wave, daughters, x̂, fftPlan, ::Union{Father, Dirac})
analyticTransformComplex!(wave, daughters, x̂, fftPlan, ::NoAve)

3. I don't find any usages for father defined in createWavelets.jl, and there is an empty averagingTypes.jl file. Is this part not complete, or simply now at

https://github.com/UCD4IDS/ContinuousWavelets.jl/blob/2ce86e57000d5b3ae70fa5a3697b25556e754d37/src/waveletTypes.jl#L129

@dsweber2 @wchak

[dsweber2]

1. Is averaging already done when doing CWT if averagingType is either Father or Dirac?

You may have noticed that daughters[:,1] is treated differently in analyticTransformComplex! https://github.com/UCD4IDS/ContinuousWavelets.jl/blob/9aceade7808f3b8748cc00c86679680f51192496/src/apply.jl#L138 That is because it is the father wavelet. In the case of analytic wavelets, this package uses the fact that most of the wavelets only act on the positive frequencies to increase performance, but the father wavelet is a real valued function, so this trick doesn't apply.

2. It seems like averagingType is simply used for multiple-dispatch in functions like analyticTrasnformReal!. In this case, maybe we can just keep the type in the arguments, something like

Yes that would be a reasonable formatting change. I wasn't aware of that syntax.

3. I don't find any usages for father defined in createWavelets.jl, and there is an empty averagingTypes.jl file. Is this part not complete, or simply now at

it definitely is? https://github.com/UCD4IDS/ContinuousWavelets.jl/blob/9aceade7808f3b8748cc00c86679680f51192496/src/createWavelets.jl#L222

I should probably delete averagingTypes.jl, as its a hold-over from a refactor that ended up being too short to be worthwhile as a standalone file (the averagingTypes are only lines 129-144 of waveletTypes.jl after all) https://github.com/UCD4IDS/ContinuousWavelets.jl/blob/2ce86e57000d5b3ae70fa5a3697b25556e754d37/src/waveletTypes.jl#L130

[dsweber2]

I just pushed a version of the crossSpectrum and waveletCoherence. There are some issues with waveletCoherence if the cwt is approximately zero, but that's to be expected.

[henry2004y]

One thing I'm not sure is the dimension of the output:

using Wavelets, ContinuousWavelets

n = 128

X = testfunction(n, "Doppler")
Y = testfunction(n, "Doppler")

c = wavelet(Morlet(π), β=2)

waveletCoherence(X, Y, c) # 128*17*1*1

For example, if the two input time-series are simply vectors, the coherence output is given as 4D, with the 3rd and 4th dimension being singleton. I briefly scanned through the code, and (correct me if I'm wrong) the current implementation allows multi-trial time-series

Z = [X Y] # pretend to have 2 trials
waveletCoherence(Z, Z, c) # 128*17*2*2

then the 3rd and 4th dimension corresponds to the coherence across different pairs of trials?

I even see condition like ndims(X) > 2 in sharedCrossSpectrum, so it supports even higher dimensional data? What does it mean?

Well, if I think in the general case this makes sense, but with simply two vectors the 4D output feels less intuitive. I would suggest separate the vector case from the general array case, such that the simplest output is a 2D time-frequency array of coherences.

---

P.S. Probably we need to also add tests for this. Am I asking for too much 🥲

[dsweber2]

3rd and 4th dimension corresponds to the coherence across different pairs of trials?

yep

I even see condition like ndims(X) > 2 in sharedCrossSpectrum, so it supports even higher dimensional data?

If a user has their data stored with a lot of extra meta dimensions, this is built to handle that (It was a use case I needed for cwt at one point). For waveletCoherence because there are two collections, I only added support for n×m, because the order of extra dimensions is pretty non-trivial. It might be better to have it throw an error for ndims(X)>2 for sharedCrossSpectrum now that I think about it.

I would suggest separate the vector case from the general array case, such that the simplest output is a 2D time-frequency array of coherences.

Yeah, I should probably include a dropdims in the output. I think there's one of those in cwt

Probably we need to also add tests for this.

I added some, such as checking that the wavelet coherence of random noise against itself is identically 1.

Am I asking for too much

I'd ignore it if I thought so. Its useful to have feedback from someone new using the package.

[henry2004y]

While testing for waveletCoherence, I get values larger than 1:

using PyPlot, Wavelets, ContinuousWavelets

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

##########

n = 2047

t = range(0, n/1000, length=n) # 1kHz sampling rate

X = testfunction(n, "Doppler")
Y = X .+ 0.2
#Y = circshift(X, 10)

c = wavelet(Morlet(π), β=2)

wc = waveletCoherence(X, Y, c) # 128*17*1*1

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

ax[1].plot(t, X)
ax[1].plot(t, Y)
ax[1].set_title("Doppler")

c = wavelet(Morlet(π), β=2)
res = ContinuousWavelets.cwt(X, 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))

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, wc[:,:,1,1]', shading="nearest")

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

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

---

Do I misunderstand coherence or is there a bug?

[dsweber2]

Because we're dividing by something that may be near zero, its a numerically unstable procedure; this is going to be particularly bad for artificial signals which are intentionally zero for a large frequency range. I implemented a test on random input, and it works fine.