为什么scipy的gaussian_filter1d无法正确计算导数?

lb3vh1jj  于 2022-11-09  发布在  其他
关注(0)|答案(1)|浏览(198)

我尝试使用关键字order来计算一个函数的导数,但是结果不正确。相反,如果我首先对函数应用高斯滤波器,然后用有限差分对它进行微分,它就可以工作了。
为了清楚起见,我想求两次微分的函数是 * 位置 *,用它我可以得到 * 加速度 *。

代码:

import numpy as np
import matplotlib.pyplot as plt
from scipy import signal
from scipy.ndimage import gaussian_filter1d

# Initial acceleration

rng = np.random.default_rng(1)
acc = np.cumsum(rng.standard_normal(2000))

# Calculate position

pos = np.zeros(len(acc)+2)
for t in range(2,len(acc)):
    pos[t] = 2*pos[t-1] - pos[t-2] + acc[t-2]

# Gaussian windows

sigma = 2
truncate = 5

acc_gaus = gaussian_filter1d(acc, sigma = sigma, truncate = truncate, order = 0)

pos_gauss = gaussian_filter1d(pos, sigma = sigma, truncate = truncate, order = 0)
acc_new_dif = pos_gauss[2:] - 2*pos_gauss[1:-1] + pos_gauss[:-2]

acc_new_gaus = gaussian_filter1d(pos, sigma = sigma, truncate = truncate, order = 2)[1:-1]

# Values

plt.figure()
plt.plot(acc_gaus[:-100], label = 'Real acceleration', alpha = 0.5)
plt.plot(acc_new_gaus[:-100], label = 'Gaussian window 2nd order', alpha = 0.5)
plt.plot(acc_new_dif[:-100], label = 'Finite differences 2nd order', alpha = 0.5)
plt.legend()
plt.ylabel('Acceleration')
plt.xlabel('Time')
plt.savefig('fig1.png', dpi = 300)

# Errors

plt.figure()
plt.plot((acc_gaus - acc_new_gaus)[:-100], label = 'Error gaussian window 2nd order')
plt.plot((acc_gaus - acc_new_dif)[:-100], label = 'Error finite differences 2nd order')
plt.legend()
plt.ylabel('Error')
plt.xlabel('Time')
plt.savefig('fig2.png', dpi = 300)

输出:
**问题:**为什么不能正常工作?在哪些情况下scipy.ndimage.gaussian_filter1d无法计算导数?

可能与以下因素有关:Does gaussian_filter1d not work well in higher orders?

i2loujxw

i2loujxw1#

这是因为高斯核的大小与输入的大小不同,如果您想获得更一致的结果,可以增加截断值,它会更接近您预期的结果。误差是累积的。使用截断值10,您可能会得到如下结果:

相关问题