In rotating machinery analysis, tachometer signals are often collected. If a 1 pulse per revolution tachometer is sampled at the same frequency as the vibration signal, then there is some jitter in the tachometer signal. Smoothing cubic splines are a preferred method of removing this jitter. This makes the data being smoothed one dimensional but not uniformly spaced.
I notice that Julia has a number of packages that will perform smoothing splines. Interpolations.jl, Splines.jl, Dierckx.jl, jipolanco/Bsplinekit.jl, SmoothingSplines, Cubicsplines, and maybe more. I would prefer a pure Julia solution. I thought someone might be able to guide me as to which one/ones I should start with.
From recent posts here in discourse, the BSplineKit.jl package seems to be very capable, it is a pure Julia solution and its author @jipolanco is actively supporting it. I would recommend starting there.
Unfortunately, smoothing splines are not currently implemented in BSplineKit.jl – only non-smoothing interpolation is included for now. I’m planning to add this in the future.
If you are looking for a pure Julia implementation, SmoothingSplines.jl seems to be the best option right now.
Couldn’t you just do a least-squares fitting using basis functions evaluated on a grid? From the docs of BSplineKit, it is my impression that all the ingredients are available.
Yes, that should be very easy to do. The only question is how to choose the spline knots, which should be sparser than the input data. Perhaps one solution would be to let the user choose the knots (as it seems to be done in Dierckx.jl).
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.
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")
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.
@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.
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.
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.
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. 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.
Also note that the front and back drums are nominally the same diameter, thus nominally the same speed.