Complex morlet wavelet (cmor) bugged ?

[Pierre-Bartet]

When using cwt on a simple sinusoid, cmor gives a weird result :

While non complex morlet or whatever other complex ones such as cgaus gives somethings that seems more correct :

[grlee77]

Can you provide a minimal Python code to reproduce this?

[Pierre-Bartet]

import numpy as np
import matplotlib.pyplot as plt
import pywt

T1 = 50
t = np.arange(0, 20*T1)
sig = np.sin(2*np.pi*t/T1)

fig, ax = plt.subplots(figsize=(12,6))
ax.plot(t, sig)

[cfs,freqs] = pywt.cwt(sig,scales=np.arange(1, 128), wavelet='cmor')
amp = np.abs(cfs)

fig, ax = plt.subplots(figsize=(12,12))
ax.imshow(amp, interpolation='nearest', aspect='auto' )

[cfs,freqs] = pywt.cwt(sig,scales=np.arange(1, 128), wavelet='morl')
amp = np.abs(cfs)

fig, ax = plt.subplots(figsize=(12,12))
ax.imshow(amp, interpolation='nearest', aspect='auto' )

[grlee77]

I agree that this does not look right.
@holgern, do you have any idea what might be going on here?

looking into it a bit more, I am not sure why cmor can be specified without any indication of what the complex Morlet parameters are. For instance, in Matlab's toolbox 'cmor' alone is not a valid wavelet name. One would have to specify a specific case such as cmor1.5-1 (although that does not currently work here).

[grlee77]

Okay, I think that maybe it is just that the frequency parameters need to be set to a different value than the default. It seems that cmor, shan and fbsp Wavelet objects have both a center_frequency and bandwidth_frequency attribute. So the equivalent of Matlab's cmor1.5-1 in PyWavelets would be:

w = pywt.ContinuousWavelet('cmor')
w.bandwidth_frequency=1.5
w.center_frequency=1

then the transform

[cfs,freqs] = pywt.cwt(sig,scales=np.arange(1, 128), wavelet=w)
fig, ax = plt.subplots(figsize=(12,12))
ax.imshow(cfs.real, interpolation='nearest', aspect='auto' )

[grlee77]

I don't think there is a bug in the transform as I see the same kind of spectrum when comparing to Matlab's cwt. To get a result in Matlab that matches the default cmor in PyWavelets, one has to set 'cmor1.0-0.5'.

I think we should modify the PyWavelets ContinuousWavelet object code to take string input such as cmor1.5-1 and extract the bandwidth/center frequency parameters for better Matlab compatibility.

Certainly we could also use more documentation and examples of the usage conventions for the CWT routines.

[Pierre-Bartet]

Thanks !
By the way it is not clear for me what is wrong with the default values in this case and why the cmor gives such a weird result compare to mor despite having the same parameters.
Also going from default values :

w.bandwidth_frequency=1.0
#w.center_frequency=0.5

to

w.bandwidth_frequency=2.0
w.center_frequency=1.0

is just a factor of 2, so I would expect the same results with a factor of 2 on the scale. In fact that is what we see on the plots : largest coefficients are at scales ~25 with default values and at ~50 otherwise.
So is it something about the numerical stability of the algorithm ? Why is only the upper triangle of the coefficient matrix affected ?

[grlee77]

So is it something about the numerical stability of the algorithm ? Why is only the upper triangle of the coefficient matrix affected ?

I am not a user of the continuous transforms and have not looked into them in detail, so I do not know why this is offhand. I can confirm that for the same sort of pattern appears in the Matlab (R2012a) output, so this behaviour is not unique to PyWavelets. Matlab result showing the magnitude of the coefficients is in the bottom panel below:

[holgern]

I found the cause of the error. It can be easily fixed.
The precision of pywt.integrate_wavelet must be increased
_cwt.py: lines 77-83

        precision = 10
        int_psi, x = integrate_wavelet(wavelet, precision=precision)
        step = x[1] - x[0]
        j = np.floor(
            np.arange(scales[i] * (x[-1] - x[0]) + 1) / (scales[i] * step))
        if np.max(j) >= np.size(int_psi):
            j = np.delete(j, np.where((j >= np.size(int_psi)))[0])

int_psi has a size of 1024. Thus for higher scales, np.max(j) >= np.size(int_psi) is true.

To fix this, precision has to be increased until np.max(j) >= np.size(int_psi) i sfalse.

Should i make a pull request, in which i try to fix this?

Should precision increased until it fits, or should set precision as function parameter?

[holgern]

    precision = 10
    int_psi, x = integrate_wavelet(wavelet, precision=precision)
    step = x[1] - x[0]
    j = np.floor(
        np.arange(scales[i] * (x[-1] - x[0]) + 1) / (scales[i] * step))
   
    if np.max(j) >= np.size(int_psi):
       precision = np.ceil(np.log2(np.max(j)))
       if precision > 25:
            precision = 25
       int_psi, x = integrate_wavelet(wavelet, precision=precision)

   if np.max(j) >= np.size(int_psi):
        j = np.delete(j, np.where((j >= np.size(int_psi)))[0])

[grlee77]

@holgern: a PR would be great!

I think accuracy is more important than speed, so doing some sort of check by default is probably best. One option would be to provide a function parameter precision that defaults to None where None means "autotune" as above. The user could then pass in a specific precision value to improve performance in cases where it is needed. Do you think that would be a reasonable choice?

[Pierre-Bartet]

It seems nice to me!

[holgern]

I investigated the problem in more detail.
It cannot be easily fixed, higher precisions does not solve the problem.

At higher scales the convolution shows artifacts:
In the following, the absolute values are plotted.
The scaled wavelet psi is at scale 100
int_psi[j.astype(np.int)][::-1]:

The coef are calculated by
coef = - np.sqrt(scales[i]) * np.diff(
np.convolve(data, int_psi[j.astype(np.int)][::-1],mode='full'))

Then the middle part is taken:

So it is not a bug, but has it cause in the convolution. The imaginary part is fine, by the way:
ax.imshow(np.imag(out), interpolation='nearest', aspect='auto' )

[holgern]

If the convolution is used with mode='valid' instead of full, the result is the following:

[holgern]

We can add the mode parameter to cwt.
mode = 'valid' means then, that only the valid part is shown and the rest is zero.
mode = 'full' means that it is tried to fill everything