0
$\begingroup$

In ssqueezepy, by finding a certain scale is it possible to do an icwt only for the scales above (or bellow) ? For exemple if I have a double sin (with $f_1$ and $f_2$), by taking the frequency $f_c$ such as $f_1$ < $f_c$ < $f_2$, going back to scales ($a_1$, $a_2$ and $a_c$), can we integrate the scales from $a_c$ (or to $a_c$) ? Would the resulting reconstructed signal would corresponds to the low frequency (or high frequency) part of my double sin ?

I tried by setting each of the elements that I didnt want to keep to zeros but the fact that it wasn't a linear scale anymore stopped it. It seems that I also can't just cut scales array because the size between scales and coefficients would be different...

What's the way to do this ?

$\endgroup$

1 Answer 1

1
$\begingroup$

I don't know if my question was unclear or trivial (or both !) but I eventually end up doing what I wanted some time ago :

So I start by creating a double sinus with two frequencies f1 and f2

plt.style.use("ggplot")

fig = plt.figure()
ax = plt.axes()

sc = 80

N1 = 2000

f1 = 1/2
f2 = 1/10

n_T = sc * f1
t_max = n_T/f1 # change here max or min if you want highest or lowest frequency to be represented on n_T periods

t1 = np.linspace(0, t_max, N1)
y1 = np.sin(2*np.pi*f1*t1) + np.sin(2*np.pi*f2*t1)

plt.plot(t1,y1, color = "firebrick", label = "test N = 2000")

plt.legend(loc = "best", frameon=True, fancybox = True,
           shadow = True, facecolor = "white")

plt.axis([-0.05*t_max, 1.05*t_max, -2, 2])

plt.title("Sine Waves")
plt.xlabel("Time")
plt.ylabel("Amplitude")

enter image description here

Then I do a cwt of that signal and represent it in frequencies

wavelet = ('morlet', {'mu': 10})

Wx2, scales2, *_ = cwt(y2, wavelet, fs=N2, scales='log')
freq2 = scale_to_freq(scales2, wavelet, N2, fs=N2/sc)
power2 = (abs(Wx2)) ** 2

plt.axhline(y=242, color='white', linestyle='--', alpha = 0.9) #for N = 2000
plt.axhline(y=168, color='white', linestyle='--', alpha = 0.9) #for N = 2000
plt.axhline(y=200, color='red', linestyle='--', alpha = 0.9)
imshow(Wx2, abs=1, yticks=freq2, title="sin with N= {}".format(N2), show=1)
print(freq2[166])
print(freq2[240])

enter image description here

Here I identify my 2 frequencies f1 and f2 and select a frequency between the two, fc (in red, $f1 > fc >f2$). I want then to split my signal in two parts, everything below fc and everything above fc. To do so I perform an icwt for the two parts selected. Here it's simply a sin for each but it would be possible to extract a low frequency in a real signal for exemple if my parts are not too close or does not intertwine frequency whise.

highfreqparticwt = icwt(Wx2[:190], wavelet, scales1[:190])
lowfreqparticwt = icwt(Wx2[190:], wavelet, scales1[190:])
plt.plot(y1, color = "blue", label = "double sin")
plt.plot(2*lowfreqparticwt, color = "red", label = "low freq part")
plt.plot(highfreqparticwt, color = "green", label = "high freq part")
plt.legend(loc = "best", frameon=True, fancybox = True,
           shadow = True, facecolor = "white")

enter image description here

Where I obtain the two parts that I wanted. I still have to multply by a constant my low frequency part to match the envelope of the original signal.

$\endgroup$

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.