I tried to use smoothness priors method from neurokit2, to detrend my signal.
https://neurokit2.readthedocs.io/en/master/_modules/neurokit2/signal/signal_detrend.html
The method was based on the research paper Method by Tarvainen et al., 2002. (Tarvainen, M. P., Ranta-Aho, P. O., & Karjalainen, P. A. (2002). An advanced detrending method with application to HRV analysis. IEEE Transactions on Biomedical Engineering, 49(2), 172-175.)
I also took a look in the paper, and the implemented code looks OK. So I really need some help on this issue:
def signal_detrend_smoothness_priors(signal, regularization=500):
N = len(signal)
identity = np.eye(N)
B = np.dot(np.ones((N - 2, 1)), np.array([[1, -2, 1]]))
import scipy as sp
D_2 = sp.sparse.dia_matrix((B.T, [0, 1, 2]), shape=(N - 2, N))
inv = np.linalg.inv(identity + regularization ** 2 * D_2.T @ D_2)
z_stat = ((identity - inv)) @ signal
trend = np.squeeze(np.asarray(signal - z_stat))
# detrend
detrended = np.array(signal) - trend
return detrended
The method constantly returns fixed two last parameters, and the signal tapers off towards the end.
End of the original signal: 1107.0, 1068.25, 1029.5, 990.75, 952.0, 939.0, 939.0, 893.0, 903.0, 1004.0, 1172.0, 1101.0, 1040.0, 986.0, 921.0, 963.0, 978.0]
End of the identified trend through the method: 192.12546294, 173.40315833, 155.18668501, 137.54799278, 120.56268372, 104.30992768, 88.87237529, 74.33606792, 60.79044746, 48.32841445, 37.04619828, 27.04344702, 18.42367655, 11.29498257, 5.76979111, 1.964643 , 963. , 978. ])
End of the detrended signal through the method: 8.32209553e+02, 8.54671586e+02, 9.66953802e+02, 1.14495655e+03, 1.08257632e+03, 1.02870502e+03, 9.80230209e+02, 9.19035357e+02, 0.00000000e+00, 0.00000000e+00])
p.s. I also tested some other code and sometimes these solutions also have issues with the last two elements in the detrended signal. -2.83746979e+01, -2.66270824e+00, 1.15376970e+02, 3.01091508e+02, 2.47593345e+02, 2.02817006e+02, 1.62618182e+02, 1.07901742e+02, -8.07242640e+02, -8.06423102e+02])
The paper says a second-order difference matrix (N-3)x(N-1). To do this you have to do:
B = np.dot(np.ones((N, 1)), np.array([[1, -2, 1]]))
D_2 = sp.sparse.dia_matrix((B.T, [0, 1, 2]), shape=(N - 3, N-1))
Now the D_2 matrix is ok (N=10, for example):
array([[ 1., -2., 1., 0., 0., 0., 0., 0., 0.],
[ 0., 1., -2., 1., 0., 0., 0., 0., 0.],
[ 0., 0., 1., -2., 1., 0., 0., 0., 0.],
[ 0., 0., 0., 1., -2., 1., 0., 0., 0.],
[ 0., 0., 0., 0., 1., -2., 1., 0., 0.],
[ 0., 0., 0., 0., 0., 1., -2., 1., 0.],
[ 0., 0., 0., 0., 0., 0., 1., -2., 1.]])