Smoothing Spline — Control Systems/Curve Fitting

Smoothing Spline

Control Systems / Curve Fitting

Smooths the last W samples with MATLAB's spaps: among all curves whose m-th derivative is square integrable, it takes the one that minimizes

ρ·E(f) + F(D^mf), with E(f) = Σ_i w_i·(y_i − f(x_i))² the weighted error at the samples and F(g) = ∫g²dx the roughness.

The weight ρ is not given – it is searched for, on every sample, until the error measure E meets the configured Tolerance. That search is what this block adds over Cubic Spline's smoothing type, which is csaps with the smoothing parameter given. The newest sample sits at x = 0 and the sample i steps back at x = −i·h, the abscissa Polynomial Fit, Fourier Fit and Spline Fit use.

Ports

• u – the sampled signal being smoothed. Scalar: one channel and its own window – see Notes.
• f – the smoothed values at the W sites of the window, [W, 1], oldest first: row W−1 is the curve at the newest sample. These are spaps's second output, and they determine the curve – it is the natural spline of order 2m through them.
• rho – the smoothing weight the search settled on, [1, 1]; spaps's third output. Zero means the error measure was already inside the tolerance and the answer is the plain weighted least-squares polynomial of degree m−1.

Parameters

• Window Length – W, how many samples the curve is fitted to. A whole number from 3 to 41, and at least m + 1. Bounded above because the two banks are unrolled at export.
• Smoothing Order – m, which derivative the roughness measures, and therefore which spline this is:
  • 1 (piecewise linear) – roughness is ∫(f′)²; the curve is a broken line through the smoothed values.
  • 2 (cubic) – ∫(f″)², the usual cubic smoothing spline, and MATLAB's default.
  • 3 (quintic) – ∫(f″′)², the quintic smoothing spline.
• Tolerance – spaps's tol, and it carries MATLAB's sign convention: a positive value is the error budget the search must meet, so bigger means smoother; a negative value skips the search and uses ρ = |tol| directly, so bigger means less smoothing. Zero is refused: in spaps it asks for the interpolating variational spline, whose answer is the data itself and whose ρ is infinite. The units are those of E – weights times the square of the signal.
• Abscissa Step – h, the spacing between consecutive samples on the fitted axis. A single positive number, and it is not cosmetic: roughness carries h^1−2m and the error measure carries h, so ρ moves roughly as h^−2m.
• Error Weights – which w_i the error measure uses:
  • Trapezoidal – MATLAB's own default when no weights are passed to spaps: the composite trapezoidal rule over the sites, so the two end samples weigh half as much as the interior ones.
  • Custom – the Custom Weights vector below.
• Custom Weights – one positive weight per sample, oldest first, or a single positive number applied to every sample (so 1 gives the unweighted error measure). Read only in the Custom mode.
• Sampling Time (s) – zero or less inherits the solver's rate; a positive value runs the block at that period.

Code export

All ten targets: Python, MATLAB, Java, Rust, C, C++, VHDL, Verilog, SystemVerilog and PLC Structured Text.

No linear system is solved in any generated core. The matrices of spaps's normal equations depend only on the sites, so the pair is diagonalized once at configuration load and each target carries two inlined banks of numbers – one that takes the window into the eigenbasis and one that brings the residual back – plus the eigenvalues. A trial value of ρ then costs one division per mode. Both banks are structural: they follow from all six settings, so re-export after changing any of them.

The search loop is bounded at 30 trial values in every target, so all ten run the same number of steps. spaps itself has no bound and was measured to need at most nine.

The three HDL targets are simulation-only real arithmetic, not synthesizable Q16.16: the search divides by a running quantity and takes a square root of one, and its trip count depends on the data. They quantize at the port boundary and in the shift register only. Note that rho leaves through a Q16.16 port there, so a configuration whose ρ runs past ±32768 – a small Abscissa Step with m = 3 reaches that easily – cannot carry it in hardware, while the seven software targets are unaffected.

Simulink bridge

None (Support::None). spaps is a Curve Fitting Toolbox function, and that toolbox ships no Simulink library at all, so there is no path a diagram could name. The bridge reports this block rather than dropping it silently, and it therefore has no parity testbench; code export verification covers it across all ten languages.

Notes

• Stateful, and discrete by nature (setDiscreteOnlyBlock(true)): the window advances once per sample. Nothing is carried between samples except the window – each window is smoothed from scratch, exactly as a call to spaps would be.
• ⚠ It is not Cubic Spline's smoothing type. That block is csaps with the smoothing parameter p given, over a table fixed in configuration, and it evaluates the curve at its input. This one searches for the smoothing weight so the error measure meets a tolerance, over a running window, at three smoothing orders, with weights – and it emits the smoothed values rather than an evaluation. Reach for Cubic Spline when the curve is known and you want its value at a point; reach for this one when the amount of smoothing should follow the data.
• ⚠ The stopping rule is MATLAB's 1 % band, so the answer is a step function of the data. spaps accepts the first trial whose error measure is within one percent of the tolerance, so two nearly identical windows can stop one iteration apart and differ by about half a percent of the residual. That is a property of the algorithm, not of any one target – but it is why a hardware target, which quantizes the input first, can disagree with the software ones on an occasional sample by more than rounding.
• The window is zero-prefilled, and the zeros count. The first W−1 outputs of a run are a startup transient, matching Polynomial Fit, Spline Fit, Detrend and Savitzky-Golay Filter. A window of zeros meets any positive tolerance, so those samples come out at ρ = 0 and f = 0.
• Scalar only. One channel and its own history; wire one block per channel.
• The fit is a compromise, and both of its ends are reachable. A tolerance larger than the error of the best degree-m−1 polynomial gives that polynomial and ρ = 0; a tolerance approaching zero drives ρ up and the curve towards interpolation.
• ⚠ Long windows at order 3 lose digits, and MATLAB loses them too. The error form's condition number grows as W^2m: at W = 13 the smoothed values reproduce spaps to 1e−13, while at the allowed maximum W = 41 with m = 3 both sides are down to about 1e−9 of each other. That is the price of the normal equations this whole family pays in exchange for putting no factorization in a generated core.
• No state space. Nonlinear in the window – the smoothing weight is found from the data – so model reduction correctly declines it.