0

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)
OverLordGoldDragon
  • 8,912
  • 5
  • 23
  • 74
Jokerp
  • 179
  • 7

1 Answers1

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

OverLordGoldDragon
  • 8,912
  • 5
  • 23
  • 74
  • 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? – Jokerp Feb 02 '23 at 14:08
  • 1
    @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. – OverLordGoldDragon Feb 02 '23 at 14:50
  • 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? – Jokerp Feb 02 '23 at 15:03
  • 1
    @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". – OverLordGoldDragon Feb 02 '23 at 15:08
  • Thank you so much for your responses. I see! So what you say is that there is two things to consider. In case i want the fluctuations as function of time for a given scale I could just Wx[scale_idx].real, however, 1) When using CWT the fluctuations would beimperfectly localized in time 2) I would need to normalize by a costant to get the desired result? What you also suggest is that issq_cwt would amend for the first issue but not the second one? Would you mind providin a MRE to do what you suggested using issq_cwt? – Jokerp Feb 02 '23 at 15:49
  • 1
    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). – OverLordGoldDragon Feb 02 '23 at 16:07