Scalograms in python

Question

I am reading this paper to learning basic concepts of dsp and I want to reproduce the following scalogram of a test signal (fig 4.2 of the paper):

It has been produced from the discretization of the formula: $$ \frac{1}{s} \int_{0}^{\Omega} f(t) \overline{g\left(\frac{t-\tau}{s}\right)} d t $$

Where $g(t)$ is the wavelet:$$ g(t)=w^{-1} e^{-\pi(t / w)^{2}} e^{i 2 \pi \eta t / w} $$

And $f(t)$ is the function: $$ \begin{array}{l} \sin \left(2 \pi \nu_{1} t\right) e^{-\pi[(t-0.2) / 0.1]^{10}} \\ \quad+\left[\sin \left(2 \pi \nu_{1} t\right)+2 \cos \left(2 \pi \nu_{2} t\right)\right] e^{-\pi[(t-0.5) / 0.1]^{10}} \\ \quad+\left[2 \sin \left(2 \pi \nu_{2} t\right)-\cos \left(2 \pi \nu_{3} t\right)\right] e^{-\pi[(t-0.8) / 0.1]^{10}} \end{array} $$

I have tried to reproduce them using the function scipy.signal.ctw(). I have coded $g(t)$ as a function of t and $\omega $. Nevertheless, I get only non-zero coefficients for the last part of the wave, I guess that I could have something wrong. How can I compute the scalogram using scipy.signal.ctw()?

Answer 1

scipy's cwt is primitive and error prone; below is via ssqueezepy.cwt:

---

Code:

Note that if you seek to code the wavelet yourself, you'll need to take it to the frequency domain first (preferably analytically via Fourier transform and then sampling it, rather than via FFT), then pass it like wavelet = Wavelet(my_func); cwt(x, wavelet).

import numpy as np
from ssqueezepy import cwt
from ssqueezepy.visuals import plot, imshow

#%%# Helper fn + params #####################################################
def exp_am(t, offset):
    return np.exp(-pi*((t - offset) / .1)**10)

pi = np.pi
v1, v2, v3 = 64, 128, 32

#%%# Make `x` & plot #########################################################
t = np.linspace(0, 1, 2048, 1)
x = (np.sin(2*pi * v1 * t) * exp_am(t, .2) +
     (np.sin(2*pi * v1 * t) + 2*np.cos(2*pi * v2 * t)) * exp_am(t, .5)  + 
     (2*np.sin(2*pi * v2 * t) - np.cos(2*pi * v3 * t)) * exp_am(t, .8))
plot(x, title="x(t) | t=[0, ..., 1], %s samples" % len(x), show=1)

#%%# Take CWT & plot #########################################################
Wx, scales = cwt(x, 'morlet')
imshow(Wx, yticks=scales, abs=1,
       title="abs(CWT) | Morlet wavelet",
       ylabel="scales", xlabel="samples")