fix the phase of CWT when using complex mother wavelets

[pavleb]

When performing CWT with complex wavelets, such as cmor the wavelet coefficients were conjugated. The equation for CWT includes conplex conjugate of the mother wavelet.

A simple example can show the problem.
Let us have two sine waves simulating a simple linear system:
Uv = 0.3np.sin(2np.pi80t)
Iv = 5np.sin(2np.pi80t + np.pi/6)

The theoretical impedance of system with such input and output signals is
Z=Uv/Iv = 0.3/5.0 * exp(-j pi/6)

When using the master branch the result is just a complex conjugate thus leading to the bug. A simple test is below:

import pywt
import numpy as np

fb = 1.0
fc = 1.0
dec_factor = 1
num_cycle = 1
fmin = 70
fmax = 90

sampling_frequency = 5000.0
t = np.arange(0,1,1.0/sampling_frequency)
Uv = 0.3*np.sin(2*np.pi*80*t)
Iv = 5*np.sin(2*np.pi*80*t + np.pi/6)

dt = 1.0/sampling_frequency
wname = 'cmor{}-{}'.format(fb,fc)

scales = np.divide(sampling_frequency,np.linspace(fmax,fmin,30))
scales=scales*fc
coefI, freqsI = pywt.cwt(Iv,scales,wname,sampling_period=dt)
coefU, freqsU = pywt.cwt(Uv,scales,wname,sampling_period=dt)


ind_coi = 100 # Remove COI ad-hoc. Ignore the first and last N samples due to COI
Zv = np.divide(coefU,coefI)
print(np.mean(Zv))
Ztheory = 0.3/5.0*np.exp(-1j*np.pi/6)
print('Tehoretical Z=%f+1j%f' % (np.real(Ztheory), np.imag(Ztheory) ))

[codecov-io]

Codecov Report

Merging #439 into master will increase coverage by <.01%.
The diff coverage is 100%.

@@            Coverage Diff             @@
##           master     #439      +/-   ##
==========================================
+ Coverage   80.79%   80.79%   +<.01%     
==========================================
  Files          23       23              
  Lines        3941     3942       +1     
  Branches      458      458              
==========================================
+ Hits         3184     3185       +1     
  Misses        681      681              
  Partials       76       76

Impacted Files	Coverage Δ
pywt/_cwt.py	89.18% <100%> (+0.3%)	⬆️

---

Continue to review full report at Codecov.

Legend - Click here to learn more
Δ = absolute <relative> (impact), ø = not affected, ? = missing data
Powered by Codecov. Last update 7fdf5d1...39b6966. Read the comment docs.

[grlee77]

Thanks for reporting this. I made a tiny change to your comment to enable Python syntax hightlighting in the example you provided.

[grlee77]

We have existing tests vs. coefficients precomuted in Matlab. It did not make sense that both the original implementation and this modification should pass. I think that was an issue with the tests themselves that I have proposed a fix for in #440.

Can you just cherry-pick the commit from that PR here? If you do, I think the tests will fail.

This means, that the current behavior in PyWavelets seems to match Matlab's behavior. Do you have a reference for the desired behavior? (I am not a CWT expert). Are you sure it is not just a different convention that is being used?

[pavleb]

I've checked your tests. The issue is actually with the old Matlab's cwt function. In 2017b onwards there is a completely new cwt function. The old one had some issues. The error in a phase is very subtle. At the moment, I am looking for a way to do point by point test. In the mean while please try the Matlab code at the end of my comment (>2017b with Wavelet toolbox) that demonstrates the issue with old CWT function when using complex wavelets.

The code uses 4 different implementations of CWT

1. wt.m From Lancaster university http://www.physics.lancs.ac.uk/research/nbmphysics/diats/tfr/wt.m
2. cwtft - Fourier transform based implementation valid for a subset of complex wavelets among which is the complex Morlet wavelet
3. New CWT implementation from 2017b onwards and
4. Old one that is used by your tests.

The signals are the same sine components. Since each function has different way of padding it is not possible to do point-by-point checks. Therefore I plot the overall ration between the wavelet coefficients of both signals.

The correct value is located at 0.0520 - 0.0300i. The first 3 implementations are correct, where as the old CWT shows conjugated values.

I would suggest one more thing. Try to run your test case for cmor1.0-2.0 against the Matlab implementation.

clear
clc
fs = 5000;
t = 0:1/fs:1-1/fs;

Uv=.3*sin(2*pi*80*t);
Iv = 5*sin(2*pi*80*t+pi/6);

scales = [55.55555556 55.98455598 56.42023346 56.8627451  57.31225296 57.7689243 ...
 58.23293173 58.70445344 59.18367347 59.67078189 60.1659751  60.66945607 ...
 61.1814346  61.70212766 62.23175966 62.77056277 63.31877729 63.87665198 ... 
 64.44444444 65.02242152 65.61085973 66.21004566 66.8202765  67.44186047 ...
 68.07511737 68.72037915 69.37799043 70.04830918 70.73170732 71.42857143];
coefsU = cwt(Uv, scales, 'cmor1.0-1.0');
coefsI = cwt(Iv, scales, 'cmor1.0-1.0');


cftu = cwtft({Uv,1/fs},'scales',scales/fs,'wavelet','morl0');
cfti = cwtft({Iv,1/fs},'scales',scales/fs,'wavelet','morl0');

[wtu,fu]=cwt(Uv,fs,'amor','VoicesPerOctave',40);%,'TimeBandwidth',40);
[wti,fu]=cwt(Iv,fs,'amor','VoicesPerOctave',40);%,'TimeBandwidth',40);


[WU,freq]=wt(Uv,fs,'nv',length(scales),'fmin',70,'fmax',90,'f0',1,'Padding','none','Wavelet','Morlet');
[WI,freq]=wt(Iv,fs,'nv',length(scales),'fmin',70,'fmax',90,'f0',1,'Padding','none','Wavelet','Morlet');


fBase = 80;
[x,ind_wt] = min(abs(freq-fBase));
[x,ind_cwt_new] = min(abs(fu-fBase));
[x,ind_cwt_old] = min(abs(scal2frq(scales,'cmor1.0-1.0')*fs-fBase));
[x,ind_cwtft] = min(abs(cftu.frequencies-fBase));

close all
plot(WU(ind_wt,:)./WI(ind_wt,:),'.')
hold all
plot(wtu(ind_cwt_new,200:end-200)./wti(ind_cwt_new,200:end-200),'*')
plot(coefsU(ind_cwt_old,200:end-200)./coefsI(ind_cwt_old,200:end-200),'x')
plot(cftu.cfs(ind_cwtft,:)./cfti.cfs(ind_cwtft,:),'o')
legend({'Lancaster','New cwt','Old cwt','FFT cwt'});

[grlee77]

Thanks @pavleb, that is pretty convincing. I really appreciate you taking the time to demonstrate the issue further. I don't have a recent Matlab available at the moment to run your example, but I was aware they had revamped their CWT routines sometime around 2015.

I am in support of making this change on our end. What do you think @rgommers?

@pavleb: Can you please also add a note regarding this to doc/release/1.1.0-notes.rst ?
I would mention it both under the "Bugs Fixed" and "Backwards incompatible changes" sections. See the older files in that folder for examples of the formatting.

I think we should probably also merge #440 after updating it to take the conjugate of the values produced by Matlab's old CWT.

[rgommers]

Yes, +1 from me as well.

[pavleb]

@grlee77 I will fix the docs and let you know. I will do it in the next couple of days.

[grlee77]

Hi @pavleb, thanks for updating the docs here.

Let me know if you need help resolving the merge conflicts (CWT has been updated to allow batched 1D transforms and FFT-based CWT since the time this was originally proposed).

[pavleb]

I fixed the conflicts. Apparently some tests are failing. I will try to check the issue in the following days. Sorry for the substantial delay.

[grlee77]

pywt/pywt/tests/test_matlab_compatibility_cwt.py

Line 151 in b904642

err = coefs_pywt - coefs

By default, the tests load precomputed results from matlab R2012a, so just conjugating those results in the line above seems to work
i.e.

err = coefs_pywt - np.conj(coefs)  # conj of Matlab R2012a result to match R2017a

[grlee77]

Hi @pavleb:
I have gone ahead and added complex conjugation of the previously computed Matlab R2012a results in the test functions to account for the bug in older Matlab. Hopefully CI will pass now so we can merge this PR soon and include it in our pending 1.1 release.

In addition to the comparisons to 3rd party software done by @pavleb, I believe the changes in this PR make the implementation consistent with Eq. 2 of Torrence and Compo's classic review article on the CWT.

[grlee77]

@pavleb: can you provide your real name? (so that you can be properly credited in the 1.1.0 release notes)
You can either add a commit modifying the name map or just post it here as a comment and I will take care of it.

[pavleb]

In order to jump over the tests, I guess it is easier for you to added to the name map. My full name is Pavle Boškoski. I saw that the code has decoration utf-8 so I guess that my last name will not make any trouble :-).

[grlee77]

I tweaked the release notes a bit in 5a38a12 if anyone wants to take a quick look. If there are no further comments this should be ready to merge.

[pavleb]

I would suggest to alter the reference to Lancaster University wt.m and include a link to their github repository with the complete toolbox https://github.com/luphysics/MODA. The function wt.m is located in this repository. Everything else is OK from my side.