2

I am really confused by the function pywt.cwt, as I've not been able to get it to work. The function seems to integrate instead of differentiating. I would like to work it as the following: Example CWT, but my graph looks like this: My CWT. The idea is to integrate the raw signal (av) with cumtrapz, then differentiate with a gaussian CWT (=> S1), and then once more differentiate with gaussian CWT (=> S2).

As you can see in the pictures, the bottom peaks of the red line should line up in the valleys, but the land under the top peaks for me, and the green line should move 1/4th period to the left but moves to the right... Which makes me think it integrates for some reason.

I currently have no idea what causes this... Does anyone happen to know what is going on?

Thanks in advance!

#Get data from pandas
av = dfRange['y']

#remove gravity & turns av right way up
av = av - dfRange['y'].mean()
av = av * -1


#Filter
[b,a] = signal.butter(4, [0.9/(55.2/2), 20/(55.2/2)], 'bandpass')
av = signal.filtfilt(b,a, av)

#Integrate and differentiate av => S1
integrated_av = integrate.cumtrapz(av)
[CWT_av1, frequency1] = pywt.cwt(integrated_av, 8.8 , 'gaus1', 1/55.2)
CWT_av1 = CWT_av1[0]
CWT_av1 = CWT_av1 * 0.05

#differentiate S1 => S2
[CWT_av2, frequency2] = pywt.cwt(CWT_av1, 8.8 , 'gaus1', 1/55.2)
CWT_av2 = CWT_av2[0]
CWT_av2 = CWT_av2 * 0.8

#Find Peaks
inv_CWT_av1 = CWT_av1 * -1
av1_min, _ = signal.find_peaks(inv_CWT_av1)
av2_max, _ = signal.find_peaks(CWT_av2)

#Plot
plt.style.use('seaborn')
plt.figure(figsize=(25, 7), dpi = 300)
plt.plot_date(dfRange['recorded_naive'], av, linestyle = 'solid', marker = None, color = 'steelblue')
plt.plot_date(dfRange['recorded_naive'][:-1], CWT_av1[:], linestyle = 'solid', marker = None, color = 'red')
plt.plot(dfRange['recorded_naive'].iloc[av1_min], CWT_av1[av1_min], "ob", color = 'red')
plt.plot_date(dfRange['recorded_naive'][:-1], CWT_av2[:], linestyle = 'solid', marker = None, color = 'green')
plt.plot(dfRange['recorded_naive'].iloc[av2_max], CWT_av2[av2_max], "ob", color = 'green')
plt.gcf().autofmt_xdate()
plt.show()
Razneesh
  • 1,147
  • 3
  • 13
  • 29
NoTimeLeft
  • 29
  • 2
  • 2
    Can you post a complete example so that one can reproduce your problem without having to guess the missing parts? Imports, dataframe, etc. – Keldorn Feb 08 '22 at 12:46
  • @user9155899 It seems that you would be better off passing output='sos' to signal.butter then signal.sosfilt(sos, av) to filter. There is a warning in the docs that output 'ba' is for backward compatibility but 'sos' should be used. Note it's difficult to help since we don't have any data and the code is not complete. https://docs.scipy.org/doc/scipy/reference/generated/scipy.signal.butter.html – Jacques Gaudin Feb 10 '22 at 15:45

1 Answers1

3

I'm not sure this is your answer, but an observation from playing with pywt...

From the documentation the wavelets are basically given by the differentials of a Gaussian but there is an order dependent normalisation constant.

Plotting the differentials of a Guassian against the wavelets (extracted by putting in an impulse response) gives the following:

differentials of a gaussian and the corresponding gausN cwt

The interesting observation is that the order dependent normalisation constant sometimes seems to include a '-1'. In particular, it does for the first order gaus1.

So, my question is, could you actually have differentiation as you expect, but also multiplication by -1?

Code for the graph:

import numpy as np
import matplotlib.pyplot as plt

import pywt

dt = 0.01
t = dt * np.arange(100)

# Calculate the differentials of a gaussian by quadrature:
# start with the gaussian y = exp(-(x - x_0) ^ 2 / dt)
ctr = t[len(t) // 2]
gaus = np.exp(-np.power(t - ctr, 2)/dt)
gaus_quad = [np.gradient(gaus, dt)]
for i in range(7):
    gaus_quad.append(np.gradient(gaus_quad[-1], dt))

# Extract the wavelets using the impulse half way through the dataset
y = np.zeros(len(t))
y[len(t) // 2] = 1
gaus_cwt = list()
for i in range(1, 9):
    cwt, cwt_f = pywt.cwt(y, 10, f'gaus{i}', dt)
    gaus_cwt.append(cwt[0])

fig, axs = plt.subplots(4, 2)

for i, ax in enumerate(axs.flatten()):
    ax.plot(t, gaus_cwt[i] / np.max(np.abs(gaus_cwt[i])))
    ax.plot(t, gaus_quad[i] / np.max(np.abs(gaus_quad[i])))
    ax.set_title(f'gaus {i+1}', x=0.2, y=1.0, pad=-14)
    ax.axhline(0, c='k')
    ax.set_xticks([])
    ax.set_yticks([])
wtw
  • 678
  • 4
  • 9