2

我对函数 pywt.cwt 感到非常困惑，因为我无法让它工作。该功能似乎是整合而不是分化。我想按以下方式工作：示例 CWT，但我的图表如下所示：我的 CWT。想法是将原始信号 (av) 与 cumtrapz 积分，然后用高斯 CWT (=> S1) 进行微分，然后再用高斯 CWT (=> S2) 进行微分。

正如你在图片中看到的，红线的底峰应该在山谷中排列，但对我来说，顶峰下方的土地，绿线应该向左移动 1/4 周期但向右移动...这让我认为它出于某种原因而集成。

我目前不知道是什么原因造成的......有没有人碰巧知道发生了什么？

提前致谢！

#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()
4

1 回答 1

3

我不确定这是你的答案，但从玩pywt...

从文档来看，小波基本上是由高斯的微分给出的，但是有一个依赖于顺序的归一化常数。

绘制高斯对小波的微分（通过放入脉冲响应提取）给出以下结果：

高斯和相应的高斯N cwt的微分

有趣的观察是，依赖于顺序的归一化常数有时似乎包括一个“-1”。特别是，它适用于第一个订单gaus1。

所以，我的问题是，你真的可以像你期望的那样进行微分，还可以乘以-1吗？

图表代码：

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([])
于 2022-02-11T14:11:54.933 回答