I am trying to calculate the derivative of a function using scipy.ndimage.gaussian_filter1d
using the keyword order
but the result is not working properly. Instead if I apply first the gaussian filter to the function and then differenciate it by finite differences it works.
For ease of clarity, the function I want to differentiate two times is the position, with which I obtain the acceleration.
Code:
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)
Question: Why is this not working properly? In which situations scipy.ndimage.gaussian_filter1d
fails to calculate the derivative?
Possibly related to: Does gaussian_filter1d not work well in higher orders?
This is coming from the fact that you gaussian kernel is not the same size as your input, if you want to get more consistant result you could increase the truncate value, it gets closer to result you expect. the error is cumulative. with truncate 10 you could get results like this: