fix the phase of CWT when using complex mother wavelets (#439)

* When performing CWT the convolution should include complex conjugate of the mother wavelet. For real wavelets this has no influence.

* Added comments for the reported bugs in the release notes

* TST: take conjugate of precomputed Matlab results from R2012a

* fix typo in Matlab reference data generation script

* DOC: update release notes

add Pavle to the authors.py name map [skip ci] [ci skip]

DOC: add link to Lancaster University MODA toolbox [ci skip]

Author: pavleb

doc/release/1.1.0-notes.rst

@@ -17,9 +17,25 @@ Deprecated features
 Backwards incompatible changes
 ==============================
 
+When using complex-valued wavelets with the ``cwt``, the output will now be
+the complex conjugate of the result that was produced by PyWavelets 1.0.x.
+This was done to account for a bug described below. The magnitude of the
+``cwt`` coefficients will still match those from previous releases.
+
 Bugs Fixed
 ==========
 
+For a ``cwt`` with complex wavelets, the results in PyWavelets 1.0.x releases
+matched the output of Matlab R2012a's ``cwt``. Howveer, older Matlab releases
+like R2012a had a phase that was of opposite sign to that given in textbook
+definitions of the CWT (Eq. 2 of Torrence and Compo's review article, "A
+Practical Guide to Wavelet Analysis"). Consequently, the wavelet coefficients
+were the complex conjugates of the expected result. This was validated by
+comparing the results of a transform using ``cmor1.0-1.0`` as compared to the
+``cwt`` implementation available in Matlab R2017b as well as the function
+``wt.m`` from the Lancaster University Physics department's
+`MODA toolbox <https://github.com/luphysics/MODA>`_
+
 Other changes
 =============
 

pywt/_cwt.py

@@ -124,6 +124,7 @@ def cwt(data, scales, wavelet, sampling_period=1., method='conv', axis=-1):
     out = np.empty((np.size(scales),) + data.shape, dtype=dt_out)
     precision = 10
     int_psi, x = integrate_wavelet(wavelet, precision=precision)
+    int_psi = np.conj(int_psi) if wavelet.complex_cwt else int_psi
 
     # convert int_psi, x to the same precision as the data
     dt_psi = dt_cplx if int_psi.dtype.kind == 'c' else dt

pywt/tests/data/generate_matlab_data_cwt.py

@@ -37,7 +37,7 @@
 try:
     all_matlab_results = {}
     for wavelet in wavelets:
-        w = pywt.Wavelet(wavelet)
+        w = pywt.ContinuousWavelet(wavelet)
         if np.any((wavelet == np.array(['shan', 'cmor'])),axis=0):
             mlab.set_variable('wavelet', wavelet+str(w.bandwidth_frequency)+'-'+str(w.center_frequency))
         elif wavelet == 'fbsp':

pywt/tests/test_matlab_compatibility_cwt.py

@@ -147,6 +147,11 @@ def _check_accuracy(data, w, scales, coefs, wavelet, epsilon):
     # PyWavelets result
     coefs_pywt, freq = pywt.cwt(data, scales, w)
 
+    # coefs from Matlab are from R2012a which is missing the complex conjugate
+    # as shown in Eq. 2 of Torrence and Compo. We take the complex conjugate of
+    # the precomputed Matlab result to account for this.
+    coefs = np.conj(coefs)
+
     # calculate error measures
     err = coefs_pywt - coefs
     rms = np.real(np.sqrt(np.mean(np.conj(err) * err)))

util/authors.py

@@ -30,6 +30,7 @@
     u('Helder'): u('Helder Oliveira'),
     u('Kai'): u('Kai Wohlfahrt'),
     u('asnt'): u('Alexandre Saint'),
+    u('pavleb'): u('Pavle Boškoski'),
 }