Inverse continuous wavelet transform off by constant factor in the y axis

Question

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()

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?

Answer 1

.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.