Smoothing interpolation for tachometer data

[rafael.guerra]

See below results for the pure Julia packages:
(1) BSplineKit.jl’s nice solution given by @jipolanco
(2) SmoothingSplines.jl
(3) Loess.jl suggested by @marius311

They work relatively well for your dataset with following caveats:
(1) allows defining the number of knots but is a bit less smooth than (2)
(2) allows defining the smoothness (all input data abscissas become knots)
(3) struggles to follow the data at start and end locations

The best is to use InspectDR.jl Plots’ backend to interactively zoom and inspect the results.

x1 =  midtime1[:]   # imported from Matlab *.mat file using MAT.jl
y1 = rotationalspeed1[:]  # imported from Matlab *.mat file using MAT.jl

using SmoothingSplines
# REINSCH, C. [1967] Smoothing by Spline Functions. Numerische Mathematik 10: 177-183
# http://eudml.org/doc/131782
spl = fit(SmoothingSpline, x1, y1, 1.0) # λ=1.0
ys1 = SmoothingSplines.predict(spl,x1) # fitted vector

using BSplineKit, LinearAlgebra, SparseArrays
# Create B-spline basis of order 4 (cubic splines) with chosen breakpoints
xbreaks = LinRange(extrema(x1)..., 200)
B = BSplineBasis(BSplineOrder(4), xbreaks)
# Construct spline from the given data by @jipolanco
C = collocation_matrix(B, x1, SparseMatrixCSC{Float64})  # evaluates basis functions at data points
coefs = qr(C) \ y1       # spline coefficients
fspline = Spline(B, coefs)

using Loess
model = loess(x1, y1, span=0.02)
yl1 = Loess.predict(model, x1)

using Plots; inspectdr()
inspectdr(size=(1000,800), lw=3, markerstrokewidth=0, ms=7, legend=:top)
Plots.scatter(x1, y1, mc=:blue,label="rotationalspeed1")
Plots.plot!(x1, ys1, lc=:red,label="(Cubic) SmoothingSpline.jl")
Plots.plot!(x1, fspline.(x1); lc=:cyan,ls=:dash, label = "Cubic BSplineKit.jl by@jipolanco")
Plots.plot!(x1, yl1, ls=:dot, lc=:black,label="Loess.jl")

Whole tachometer signal:

Zoom-in at start:

Zoom-in at peak:

Zoom-in at end:

[Jake]

Thank you very much Rafael. That is extremely helpful. And it appears that you have spent some time trying to optimize the parameters for each of the smoothing algorithms. I modified your code a little to include the reading of the data file and time each of the smoothing algorithms. I found it odd that the type of the variable from the MATread had to be changed with the x1 = midtime1[:] command, but it worked nicely.

using Plots
using MAT
using BenchmarkTools
d = matread("tachdata.mat")
midtime1 = d["midtime1"]; rotationalspeed1 = d["rotationalspeed1"]
midtime2 = d["midtime2"]; rotationalspeed2 = d["rotationalspeed2"]
midtime3 = d["midtime3"]; rotationalspeed3 = d["rotationalspeed3"]
midtime4 = d["midtime4"]; rotationalspeed3 = d["rotationalspeed3"]

plot(d["midtime1"], d["rotationalspeed1"], label = "paper roll")
plot!(d["midtime2"], d["rotationalspeed2"], label = "back drum")
plot!(d["midtime3"], d["rotationalspeed3"], label = "front drum")
plot!(d["midtime4"], d["rotationalspeed4"], label = "unwind roll")

x1 =  midtime1[:]   # imported from Matlab *.mat file using MAT.jl
y1 = rotationalspeed1[:]  # imported from Matlab *.mat file using MAT.jl

using SmoothingSplines
println("SmoothingSplines")
# REINSCH, C. [1967] Smoothing by Spline Functions. Numerische Mathematik 10: 177-183
# http://eudml.org/doc/131782
@btime spl = fit(SmoothingSpline, x1, y1, 1.0) # λ=1.0
@btime ys1 = SmoothingSplines.predict(spl,x1) # fitted vector

using BSplineKit, LinearAlgebra, SparseArrays
println("BSplineKit")
# Create B-spline basis of order 4 (cubic splines) with chosen breakpoints
@btime xbreaks = LinRange(extrema(x1)..., 200)
@btime B = BSplineBasis(BSplineOrder(4), xbreaks)
# Construct spline from the given data by @jipolanco
@btime C = collocation_matrix(B, x1, SparseMatrixCSC{Float64})  # evaluates basis functions at data points
@btime coefs = qr(C) \ y1       # spline coefficients
@btime fspline = Spline(B, coefs)
@btime yb1 = fspline.(x1)

using Loess
println("Loess")
@btime model = loess(x1, y1, span=0.02)
@btime yl1 = Loess.predict(model, x1)

using Plots; inspectdr()
inspectdr(size=(1000,800), lw=3, markerstrokewidth=0, ms=7, legend=:top)
Plots.scatter(x1, y1, mc=:blue,label="rotationalspeed1")
Plots.plot!(x1, ys1, lc=:red,label="(Cubic) SmoothingSpline.jl")
Plots.plot!(x1, yb1; lc=:cyan,ls=:dash, label = "Cubic BSplineKit.jl by@jipolanco")
Plots.plot!(x1, yl1, ls=:dot, lc=:black,label="Loess.jl")

matwrite("tachplot.mat", Dict(
    "x" => x1,
     "y" => hcat(y1, ys1, yb1, yl1)))

The timing results generally give the time to calculate the spline and then the time to calculate the curve from the spline with more operations for BSplineKit and are shown below:

julia> include(“splinefittachdata.jl”)

SmoothingSplines
1.771 ms (26926 allocations: 3.12 MiB)
400.900 μs (2 allocations: 52.58 KiB)

BSplineKit
16.000 μs (4 allocations: 112 bytes)
71.473 ns (1 allocation: 64 bytes)
16.195 ms (32 allocations: 1.00 MiB)
1.903 ms (137 allocations: 6.38 MiB)
74.089 ns (1 allocation: 64 bytes)
444.100 μs (3 allocations: 52.66 KiB)

Loess
817.901 ms (7416603 allocations: 714.19 MiB)
32.767 ms (717141 allocations: 25.95 MiB)

This shows that SmoothingSplines has the fastest execution time as well as giving really good performance.

This was done on Win10 with a i7 9th generation processor.

Two issues:

When I execute the script with include(“splinefittachdata.jl”) I don’t get any plots. Is there something that I am missing?

When I plot from the REPL I was not able to zoom in on the plots so I cheated and used Matlab (Thus the writing of data to Matlab). Any suggestions on what that problem may be would be appreciated. I am working on installing on a Linux machine in case that gives a different result.

[JM_Beckers]

Another pure Julia smoothing interpolation tool (in N dimensions), mostly used in oceanography but it can work with any data if you use the core routine divandrun:

https://github.com/gher-ulg/DIVAnd.jl

[rafael.guerra]

@Jake, concerning the first issue, could you try in that case to include the plot commands inside display()?

Regarding the second issue, I have tried the REPL in Win10 Julia 1.6.0-beta1.0, and everything plots and zooms like a charm with inspectdr() backend. Please check this page for the: InspectDR: Mouse/Keybindings

[Jake]

Wrapping each plot command within display() worked, thanks.

I tried the program but replacing the first tach with the second tach. This caused BSplineKit to error with the following issue

julia> coefs = qr(C) \ y1
ERROR: DimensionMismatch("size(F) = (6717, 202) but size(B) = (5987,)")
Stacktrace:
 [1] _ldiv_basic(F::SuiteSparse.SPQR.QRSparse{Float64, Int64}, B::Vector{Float64})
   @ SuiteSparse.SPQR C:\buildbot\worker\package_win64\build\usr\share\julia\stdlib\v1.6\SuiteSparse\src\spqr.jl:393
 [2] \(F::SuiteSparse.SPQR.QRSparse{Float64, Int64}, B::Vector{Float64})
   @ SuiteSparse.SPQR C:\buildbot\worker\package_win64\build\usr\share\julia\stdlib\v1.6\SuiteSparse\src\spqr.jl:433
 [3] top-level scope
   @ REPL[207]:1

I am not sure why either.

[rafael.guerra]

For second series you need to change also the abscissas. Keeping the same variable names this means:

x1 = midtime2[:] 
y1 = rotationalspeed2[:]

[Jake]

That is what I did and SmoothingSplines and Loess work, though SmoothingSplines seems better and faster. It looks like I will work with that package.

[rafael.guerra]

I still think you may need to check your variables definition as the BSplineKit package runs fine for me on your tach data. See plot below for the second tach:

[Jake]

With some more experimentation, there is an intermittent problem when I use @btime on front of the statements. When I remove this macro the code seems to work fine. It is limited experimentation so far. And it is using 1.6.0-rc0. Do you think it could be a bug in this version?

I figured out how to zoom in the plot, I need to click and hold the right mouse button. I have not yet figured out how to unzoom or reset the zoom level. I have also not found a reference as to how to do this.

The zooming is really nice though, thanks for pointing it out.

[rafael.guerra]

As posted above, please check this page for the: InspectDR: Mouse/Keybindings