问为什么gaussian_filter1d的导数计算不正确？
EN

Stack Overflow用户

提问于 2022-11-04 09:59:17

回答 1查看 30关注 0票数 0

我试图使用scipy.ndimage.gaussian_filter1d使用关键字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)

输出：