0
$\begingroup$

I have implemented the Continuous Wavelet Transform using the pycwt library(https://github.com/regeirk/pycwt/blob/master/pycwt/wavelet.py) and its inverse using Morlet wavelets, however, upon calculating the inverse, the signal generated is accurate but with a constant discrepancy.

The code to produce icwt on pycwt library had some issues that i tried to fix, which improved the result significantly. This is the code I use for the ICWT.

import pycwt

def _check_parameter_wavelet(wavelet):
    mothers = {'morlet': pycwt.mothers.Morlet}
    # Checks if input parameter is a string. For backwards
    # compatibility with Python 2 we check either if instance is a
    # `basestring` or a `str`.
    try:
        if isinstance(wavelet, basestring):
            return mothers[wavelet]()
    except NameError:
        if isinstance(wavelet, str):
            return mothers[wavelet]()
    # Otherwise, return itself.
    return 

def icwt(W, sj, dt, dj=1/12, wavelet='morlet'):
    """Inverse continuous wavelet transform.
    Parameters
    ----------
    W : numpy.ndarray
        Wavelet transform, the result of the `cwt` function.
    sj : numpy.ndarray
        Vector of scale indices as returned by the `cwt` function.
    dt : float
        Sample spacing.
    dj : float, optional
        Spacing between discrete scales as used in the `cwt`
        function. Default value is 0.25.
    wavelet : instance of Wavelet class, or string
        Mother wavelet class. Default is Morlet
    Returns
    -------
    iW : numpy.ndarray
        Inverse wavelet transform.
    Example
    -------
    >> mother = wavelet.Morlet()
    >> wave, scales, freqs, coi, fft, fftfreqs = wavelet.cwt(var,
           0.25, 0.25, 0.5, 28, mother)
    >> iwave = wavelet.icwt(wave, scales, 0.25, 0.25, mother)
    """
    wavelet = _check_parameter_wavelet(wavelet)

    a, b = W.shape
    c = sj.size
    if a == c:
        sj = (np.ones([b, 1]) * sj).transpose()
    elif b == c:
        sj = np.ones([a, 1]) * sj
    else:
        raise ValueError('Input array dimensions do not match.')

    # As of Torrence and Compo (1998), eq. (11)
    iW = (dj * np.sqrt(dt) / (wavelet.cdelta * wavelet.psi(0)) *
          (np.real(W) / np.sqrt(sj)).sum(axis=0))

    return iW

Here I providea minimum example where I do the CWT and reconstruct with my ICWT and the original ICWT. However:

# Parameters
N  = 20000
dt = 1
dj = 1/4

#Create a signal forexample
sig = np.sin(np.linspace(0,  5*np.pi, N))+np.random.rand(N)

# DO the CWT analysis
W, sj, freqs, coi,_,_ = pycwt.cwt(sig, dt, dj, wavelet='morlet')

# original icwt
iwave                 = pycwt.icwt(W, sj, dt, dj, 'morlet')

# fixed icwt
iwave_fixed           = icwt(W, sj, dt, dj, 'morlet')

#Plot results
plt.plot(iwave, label='ICWT')
plt.plot(iwave_fixed, label='fixed ICWT')
plt.plot(sig, label='Original')

plt.legend()

enter image description here

As you can see when I compute the inverse, the resulting signal is off by some constant factor (but otherwise correct). Depending on which frequencies of wavelet I use for the transforms, the resulting signal is off by a different constant factor.

Any ideas why I am getting this shift and how I could fix it?

$\endgroup$

1 Answer 1

1
$\begingroup$

.sum(axis=0) is the one-integral inverse, and depends on the forward transform. Check against the list of conditions outlined there. Here's a working icwt.

Also CWT and iCWT are independent of the sampling rate (1/dt) unless we're normalizing for something, and that can be quite tricky.

$\endgroup$
4
  • $\begingroup$ Thanky you so much for sending me your library! it looks amazing! If I may i have two questions related to your library. 1) When doing the CWT how could I extract the frequencies corresponding to each scale? 2) Would it be possible to to the ICWT for a given scale and return the fluctuations of the timeseries corresponding to this scale? $\endgroup$
    – Jokerp
    Commented Feb 1, 2023 at 16:26
  • 1
    $\begingroup$ @Jokerp Glad you like it. 1) see here, 2) yes but this should be asked separately $\endgroup$ Commented Feb 1, 2023 at 16:56
  • $\begingroup$ Thank you so much! I will create a quesiton later today! $\endgroup$
    – Jokerp
    Commented Feb 1, 2023 at 17:01
  • $\begingroup$ I asked the question here: dsp.stackexchange.com/questions/86455/…. it would really help if you could provide your input. $\endgroup$
    – Jokerp
    Commented Feb 1, 2023 at 19:21

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.