Inverse continuous wavelet transform, how to obtain fluctuations at a given scale with ssqueezepy?

Question

How can I obtain the fluctuations of a timeseries at a specific scale using the ssqueezepy library for inverse continuous wavelet transform (ICWT)?

I have a minimum reproducible example that uses the squezepy library to perform the continuous wavelet transform (CWT) and estimate the ICWT. However, what I would like to do is to perform the ICWT for a particular scale and return the fluctuations of the timeseries corresponding to that scale. Can someone help me with this problem?

import numpy as np


import ssqueezepy
def PSD_cwt_ssqueezepy(x, y, z, dt, nv =16, scales='log-piecewise', wavelet=None, wname=None, l1_norm=False, est_PSD=False):
    """
    Method to calculate the wavelet coefficients and  power spectral density using wavelet method.
    Parameters
    ----------
    x,y,z: array-like
        the components of the field to apply wavelet tranform
    dt: float
        the sampling time of the timeseries
        
    scales: str['log', 'log-piecewise', 'linear', 'log:maximal', ...]
                / np.ndarray
            CWT scales.
    Returns
    -------
    W_x, W_y, W_zz: array-like
        component coeficients of th wavelet tranform
    freq : list
        Frequency of the corresponding psd points.
    psd : list
        Power Spectral Density of the signal.
    scales : list
        The scales at which wavelet was estimated
    """
    
    if wavelet is None:
        wavelet    = ssqueezepy.Wavelet(('morlet', {'mu': 13.4}))
    else:
        wavelet    = ssqueezepy.Wavelet((wname, {'mu': 13.4}))  
        
    if  scales is  None:
        scales = 'log-piecewise'

    # Estimate sampling frequency
    fs         = 1/dt
    # Estimate wavelet coefficients
    Wx, scales = ssqueezepy.cwt(x, wavelet, scales, fs, l1_norm=l1_norm, nv=nv)
    Wy, _      = ssqueezepy.cwt(y, wavelet, scales, fs, l1_norm=l1_norm, nv=nv)
    Wz, _      = ssqueezepy.cwt(z, wavelet, scales, fs, l1_norm=l1_norm, nv=nv)
    
    # Estimate corresponding frequencies
    freqs      = ssqueezepy.experimental.scale_to_freq(scales, wavelet, len(x), fs)
    
    if est_PSD:
        # Estimate trace powers pectral density
        PSD = (np.nanmean(np.abs(Wx)**2, axis=1) + np.nanmean(np.abs(Wy)**2, axis=1) + np.nanmean(np.abs(Wz)**2, axis=1)   )*( 2*dt)
    else:
        PSD = None
    
    return Wx, Wy, Wz, freqs, PSD, scales, wavelet


# Parameters
N       = 20000
dt      = 1
nv      =  16
scales  ='log-piecewise'
wavelet = None
wname   = None 
l1_norm = False 
est_PSD = False

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


# Do the CWT analysis to obtain wavelet coeeficients
Wx, Wy, Wz, freqs, PSD, scales,wavelet = PSD_cwt_ssqueezepy(sig, sig, sig, dt, nv , scales, wavelet, wname, l1_norm, est_PSD)

# Obtain ICWT for one of the components Wx

rec_signal = ssqueezepy.icwt(Wx, 
                  wavelet, 
                  scales,
                  nv,
                  l1_norm)

Answer 1

Per how the one-integral inverse works, this is as simple as Wx[scale_idx].real, within a scaling. However it's probably not what you want.

The CWT is inherently imperfectly localized in both time and frequency, so any result that desires said result must be handled carefully. Relevant article. If you really want contents "at" a time or frequency, you need instantaneous methods, like synchrosqueezing.

If instead what you want is icwt of an individual signal component, that's a well-defined task. Here scale_idx may vary with time, so specifying is non-trivial - refer to this post, which contains unfinished code of a tool I made to facilitate manual extraction. There's also ssqueezepy.extract_ridges.

Lastly, ssqueezepy doesn't have a built-in for inverting CWT over specific scales vs time, but we can simply use issq_cwt, and for most part the result will be within a constant of the true result. You could try combining the right parts of issq_cwt (_invert_components) and icwt to improve the results. But, if you want to invert over entire scales (not vs time), it's as simple as passing in Wx[i0:i1] and scales[i0:i1] - just careful with scales='log-piecewise', pass in only the log or linear portion at a time to avoid edge cases (use idx = logscale_transition_idx(scales) to get where the transition happens).