0
$\begingroup$

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)
$\endgroup$

1 Answer 1

1
$\begingroup$

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

$\endgroup$
6
  • $\begingroup$ Thank you for your reply. I am not exactly sure i understand your answr. What you suggest is that If I wanted to get the fluctuations at a given scale I could do this using ssq_cwt and then issq_cwt on Wx[scale_idx].real? What I am intersted in getting is several timeseries. Each corresponding to a given scale of the CWT( or ssq_CWT) that results from taking the ICWT (or issq_cwt ). Is this even possible? $\endgroup$
    – Jokerp
    Commented Feb 2, 2023 at 14:08
  • 1
    $\begingroup$ @Jokerp Wx[scale_idx].real are the fluctuations. Equals. If scale of interest is 5, and scales[20] == 5, then Wx[20].real is what you need. For several scales at once, we'd do e.g. icwt(Wx[5:30], scales[5:30]), that returns the combined signal. $\endgroup$ Commented Feb 2, 2023 at 14:50
  • $\begingroup$ Isnt Wx[scale_idx].real the fluctuation in k space? What I meant is after obtaing the wavelet coefficients ICWT them at each scale to get the fluctuation in time. Is this what you meant? $\endgroup$
    – Jokerp
    Commented Feb 2, 2023 at 15:03
  • 1
    $\begingroup$ @Jokerp The coefficient is the inverse, that's the beauty of satisfying the one-integral inverse criteria in the forward transform. It's what I describe on top here. Don't know what you mean by "k space". $\endgroup$ Commented Feb 2, 2023 at 15:08
  • 1
    $\begingroup$ That's a bit much back & forths, so briefly. Yes & yes 1 & 2. issq_cwt is for inverting ssq_cwt, not cwt, but it has _invert_components. It solves neither 1 nor 2 with cwt, but icwt solves 2 for cwt. ssq_cwt solves 1 and 2, but 1 just partially. I suggest reading through some code and experimenting with signals to better understand what's happening. -- Also I should mention that this "direct inversion" only works with l1_norm=True (default). $\endgroup$ Commented Feb 2, 2023 at 16:07

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.