0

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)

Outputs: enter image description here enter image description here

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?

Puco4
  • 491
  • 5
  • 16

1 Answers1

1

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:

enter image description here

enter image description here

amirhm
  • 1,239
  • 9
  • 12
  • So how could you know in practice whether the derivative is being calculated properly with the Gaussian filter? – Puco4 Nov 04 '22 at 13:57
  • it is calculated always correctly, what is different is that you should see what are you interested in doing precisely, when you use the order higher than 0 to use the derivatives of Gaussian and truncates means that you should obviously accumulate the error. – amirhm Nov 04 '22 at 14:12
  • But I mean how would you decide which truncate value to choose in order to keep a low accumulated error? – Puco4 Nov 04 '22 at 14:24
  • truncate value is how many times of sigma you want to keep, look at gaussian distribution you could see by truncating you are losing data https://en.wikipedia.org/wiki/Standard_error. it is discrete and finally could not be exact, it will drift, – amirhm Nov 04 '22 at 14:42
  • I don't really see why it should accumulate the error though and how to think of a value in a given situation. – Puco4 Nov 04 '22 at 15:01