# High-order derivatives of sampled data

Hi,

I have discrete data sampled at a constant rate. I need the first, second, and third derivatives with respect to time. Is there a Julia package that could help me do that? Things that I’ve considered so far:

• Use finite differences from either FiniteDiff.jl or FiniteDifferences.jl, but both appear to differentiate Julia functions, not discrete data.
• Recursively use spline interpolation and then differentiate the spline. That led me to a discontinuous second derivative, which is no good.
• Use ApproxFun.jl, but that again seems to work only for continuous functions, not discrete data.

Anyone have any better ideas?

Thanks!

Daniel

1 Like

Is O(\Delta t^2) accuracy sufficient? If so, you could just use low-order finite-difference formulas, e.g. see here.

2 Likes

Yes, I think that finite differences would work fine—I was just surprised that I couldn’t find an existing implementation in the Julia ecosystem.

Thanks!

@DanielIngraham, have you tried Dierckx.jl or BSplineKit.jl for this?

A very simple example is provided below using Dierckx, but would recommend BSplineKit as it is native Julia and greatly supported by the author:

``````using Dierckx
f(x) = exp(-2x)
x = 0:0.1:2;  y = f.(x)
spl = Spline1D(x, y; k=5)
D1f = derivative(spl, x; nu=1)
D2f = derivative(spl, x; nu=2)
D3f = derivative(spl, x; nu=3)
``````

2 Likes

You could take derivatives using a Gaussian Process regression, if you select a differentiable kernel. You could the Gaussian process julia package to do the regression.

1 Like

Sounds like using a massive sledgehammer to crack a nut, but might be wrong.

1 Like

Below is @rafael.guerra’s example using BSplineKit.

``````using BSplineKit
f(x) = exp(-2x)
x = 0:0.1:2
y = f.(x)
spl = interpolate(x, y, BSplineOrder(6))
D1f = diff(spl, Derivative(1))
D2f = diff(spl, Derivative(2))
D3f = diff(spl, Derivative(3))

using CairoMakie
D3f_exact(x) = -8f(x)
fig = Figure(resolution = (400, 300))
ax = Axis(fig[1, 1])
plot!(ax, 0..2, D3f_exact; label = "Exact")
scatter!(ax, x, D3f.(x); color = :red, label = "BSplineKit f‴(x)")
axislegend(position = :rb)
`````` Here I used B-splines of order k = 6 (piecewise polynomials of degree 5), but the order can be arbitrarily increased or decreased.

2 Likes

A simpler basis would be polynomial, which we could get using ApproxFun:

``````using ApproxFun

function vandermonde(S,n,x::AbstractVector)
V=Array{Float64, 2}(undef, length(x), n)
for k=1:n
V[:,k]=Fun(S,[zeros(k-1);1]).(x)
end
V
end

x=collect(0:0.1:2)
dom = domain(Fun(identity,minimum(x)..maximum(x)))
v=exp.(-2x)

S=Taylor()
k = 9 # polynomial order
V=vandermonde(S,k,x)
f=Fun(S,V\v)

order = 3 # derivative order
g = Derivative(dom,order) * f

plot(x,-8 .* v, label="Exact")
scatter!(x, g.(x), label="ApproxFun", legend=:bottomright)
``````

This is the least-squares fit of the data used in this thread using the Taylor polynomial. ApproxFun has a bunch of other bases. ApproxFun also includes piecewise domains for combining multiple bases, which you may be able to use if the data is piecewise continuous and you apply a change-detection algorithm.

1 Like

As an aside, regarding piecewise continuous data, e.g., jump discontinuities, I just found a really nice julia package: NoiseRobustDifferentiation. It’s based off of work from Rick Chartrand (a Los Alamos researcher):

• Rick Chartrand, “Numerical differentiation of noisy, nonsmooth data,” ISRN Applied Mathematics, Vol. 2011, Article ID 164564, 2011.

A while back I was interested in real-time numerical differentiation for nonlinear systems using embedded controllers. At the time, I developed a variant of the numerical differentiator described in:

• M. Mboup, C. Join, and M. Fliess, “Numerical differentiation with annihilators in noisy environment.” Numerical Algorithms, 50 (4), 439–467, 2009.

There’s a nice python/Matlab BSD-licensed implementation of these algebraic differentiators here. I didn’t see an implementation in julia. I may end up translating it because it allows for explicitly setting the low-pass characteristics of the differentiator at higher frequencies. This is useful if you want as output a causal LTI filter with coefficients that could implement for real-time use.

Numerical differentiation is a very interesting rabbit hole with so many ways of attacking the problem, each with its own advantages / disadvantages.

3 Likes

Excellent, thanks everyone! I ended up going with BSplineKit. This community is great.

2 Likes