Complex Morlet wavelet in CWT is not symmetric at fractional scale

[amanita-citrina]

Complex Morlet coefficients should always be symmetric

This issue may be related to issue #531.

The complex Morlet wavelet at fractional scales is not symmetric in the real part nor anti-symmetric in the imaginary part. This can be seen by computing the CWT of a Dirac signal (single impulse).

In my application, this caused strange artefacts similar to #531. They vanished when I applied the equations from Wikipedia directly. I chose the sampling instants such that the wavelet coefficients were symmetric/anti-symmetric.

To reproduce the issue, consider the code snippet below. (Integer scales yield the expected symmetries.)

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

x = np.arange(-16,16) == 0
coef,freqs = pywt.cwt(x, 2.3, 'cmor2.0-1.0')
y = np.vstack([np.real(coef), np.imag(coef)]).T
plt.plot(y)

[grlee77]

Thanks for reporting the issue. Can you perhaps share a bit more about the solution you propose? If you make a pull request (or point me to a branch where you have implemented it), even if it is not perfect, that could give a clearer start point for discussion.

Does making the one-line change to precision proposed in #531 also help alleviate the problem in your case or is it still persistent? I am guessing the issue is one of how the discretization of the continuous transform is being done.

The original contributor of CWT to this library implemented it consistently with an older release of Matlab, which appears to be suboptimal in some ways. I would be open to alternative/improved implementations, but do not personally have the free time to work on an alternative CWT implementation. Contributions in this area would be greatly appreciated.

[amanita-citrina]

Actually, I don't have a solution right now. I'll have to dig into the code and see if I can come up with a proposal. This will take a few days, though.

In general terms, I made sure that the sampling instants of the wavelet signal were symmetric by choosing an odd number of samples centered around T=0.

I guess that the current implementation has an even number of samples, and that due to some rounding or unlucky scaling the resulting sampling instants are not symmetric.

[amanita-citrina]

In _cwt.py, the wavelet is integrated (line 126):

int_psi, x = integrate_wavelet(wavelet, precision=precision)

using 2**precision samples from the wavelet. The integration simply uses np.cumsum().

Later on, a number of samples is selected for the actual convolution. As found in issue #531, the sample indices are obtained by rounding the floating point indices towards zero (floor function).

Finally, the integration is reversed in line 183:

coef = - np.sqrt(scale) * np.diff(conv, axis=-1)

A few observations:

• Increasing 'precision' seems not to help with wavelet symmetry.
• I tried to somehow 'symmetrize' the wavelet sampling instants, but without success.
• If I am not mistaken, the integrate/diff steps cancel out (up to rounding) and could be omitted.
• In this case, the wavelet should simply be evaluated at the (symmetric) scale-dependent sampling instants

Unfortunately, this is about as far as I am able to help. From the remarks in #531 I suspect that it would be a good idea to reimplement cwt() as in scipy, i.e., without integration.