High-order derivatives of sampled data

General Usage
,
1

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
2

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

2 Likes
3

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!

4

@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)

Dierckx_high_order_derivatives
Dierckx_high_order_derivatives1800×1200 146 KB

3 Likes
5

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.

https://stats.stackexchange.com/questions/373446/computing-gradients-via-gaussian-process-regression

1 Like
6

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

1 Like
7

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)

interp

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
8

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
9

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
10

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

2 Likes