Smoothing interpolation for tachometer data

Another to look at is Loess.jl, its pure Julia and always impresses me how much it puts a smooth line through the data exactly like I would have.

1 Like

Can you add a plot so we can see the noise structure of your data?

I am a bit intrigued that smoothing of tachometer data seems to have been not so often considered as statistical signal processing problem.

2 Likes

Possibly I did not understand, but please note that in Dierckx.jl, there is an option to have the number and positions of knots chosen automatically, controlled by a smoothing factor parameter s.

Yes, that’s the kind of smoothing spline that I was referring to in my first answer, in which one has a tunable smoothing parameter.

Dierckx can do both smoothing splines (with a smoothing parameter s), and what seems to be least-square fitting (the second form of Spline1D), which is what @Tamas_Papp was referring to. The latter takes a set of spline knots instead of taking a smoothing parameter.

1 Like

I would not overthink this, and just use an equispaced grid.

In that case, here is a quick example showing how such a fit would be done in the current version of BSplineKit. I may add a higher level interface for this in the future.

using BSplineKit
using LinearAlgebra
using SparseArrays
using Random
using Plots

# Exact function
f(x) = exp(-x) * cos(8x)

# Locations where data is available (here is equispaced for simplicity)
xdata = 0:0.025:1

# Spline breakpoints (should be sparser than the data points)
xbreaks = 0:0.1:1

# Data with some random noise
ydata = map(x -> f(x) + 0.1 * randn(), xdata)

# Create B-spline basis of order 4 (cubic splines) with the chosen breakpoints
B = BSplineBasis(BSplineOrder(4), xbreaks)

# Construct spline from the given data
C = collocation_matrix(B, xdata, SparseMatrixCSC{Float64})  # evaluates basis functions at data points
coefs = qr(C) \ ydata       # spline coefficients
fspline = Spline(B, coefs)

plot(f, 0, 1; label = "f(x)", xlabel = "x", linestyle = :dash)
scatter!(xdata, ydata; label = "data")
plot!(x -> fspline(x), 0, 1; label = "fit")

5 Likes

@marius311, looking at the Loess.jl page, the robust fit smooth line seems to be way off?

Tested it below against SmoothingSplines.jl. and good results were obtained using your recommended Loess span parameter.

using Loess, SmoothingSplines, Plots
xs = sort(10 .* rand(200))
ys = sin.(xs)
yn = ys .+ (rand(200) .- 0.5)
yn[96:2:104] .= 2.9;  # Add massive outliers
model = loess(xs, ys, span=0.5)
# us = range(extrema(xs)...; step = 0.1)
vs = Loess.predict(model, xs)

Plots.plot(xs, ys, label="True signal", ylims=(-2.0,3))
Plots.scatter!(xs, yn, ms=3, label="Noisy samples w/ outliers")
Plots.plot!(xs, vs, ls=:dash, label="Loess.jl")

spl = fit(SmoothingSpline, xs, ys, 1.0) # λ=1.0
yp = SmoothingSplines.predict(spl,xs) # fitted vector
Plots.plot!(xs, yp, lc=:red,ls=:dash,label="(Cubic) SmoothingSpline.jl")

Loess_vs_SmoothingSplines

2 Likes

Yea that example it looks kind of weird. The main option to play with is span, lowering it does less smoothing, with loess(xs, yn, span=0.5) your example works much better I think.

1 Like

@marius311, thanks for the feedback. It looks much better now with your recommended Loess span parameter. I have updated example above accordingly and included some massive outliers to spice it up further.

Here is a very simple constant speed example of the raw data with some spline fits that I did using Matlab and matlab central and other toolboxes.

It gives an example of the jitter in the tach signal. I will lookup some other tach signals and post them from non-constant speed sources to give some further idea.

For this curve equal knot spacing is certainly more than adequate.

1 Like

A bit more detail for those curious.

If you picture some reflective tape on a rotating shaft and then a digital optical pickup on the shaft and collect the data with a standard digital signal analyser that includes anti-alias filtering you will get a signal much like this.

.

Note the signal overshoot due to the filtering.

In my case the equipment is usually large industrial equipment that changes speed slowly. In the case of a paper machine winder, because of the fragile paper web, the acceleration rate of change (jerk) occurs over a few seconds and the speed profile looks like the following figure, where the winder speed was reduced in the middle, probably due to a paper issue. with a close up view as this.

If you look closely, you will see that the different splines tried have differing capability of following the acceleration change.

Another use case in the automotive industry would be following a car engine rotational speed as the transmission is shifting. In the ideal world, knot placement would automatically occur around the gear shift, or acceleration change, or high jerk. But the spline would smooth out the jitter in the signal. In practice, so far I have used uniformly spaced splines, though have thought about more clever knot placement.

I have done some testing where a variable speed AC drive was used for a speed sweep. In this case, it was more like a stepped speed sweep rather than a uniformly varying speed sweep. Again an ideal situation for placing knots at the locations of speed change. But again that will be much more complex than uniformly distributed knots.

1 Like

Would it be possible to put in a dropbox your tachometer signal and share it?

1 Like

This example is from a winder similar to that shown in https://www.valmet.com/board-and-paper/board-and-paper-machines/winders-and-roll-handling/winding/two-drum-winders/.

The data is available for a week at Transfer - Dropbox. Note the the data contains the raw tachometer signals which have been processed into the time and speed data which are also contained in the file and plotted using the code below.

using Plots
using MAT

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

The code and the plot should explain the variables. tachs Also note that the front and back drums are nominally the same diameter, thus nominally the same speed.

I look forward to reviewing all your comments in more detail and trying them out on this data.

@Jake, please confirm which of the datasets you would like to interpolate.

Perhaps midtime1 and rotatationalspeed1 is a good place to start.

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:

3 Likes

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.

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

2 Likes

@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

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.