2
$\begingroup$

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

enter image description here

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

$\endgroup$
0
1
$\begingroup$

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

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.