I have a lot of experience with the GCV spline code (received it in 1985 on a 9-track tape!).

The Fortran source code is still available on Netlib at

https://www.netlib.org/gcv/gcvspl. After the first page or two, look at the "SUBROUTINE GCVSPL" with good technical information. There is an accompanying paper with the full mathematical details [1]. Tests and comparisons with other methods are in [2].

GCV stands for "Generalized Cross Validation" and it is a method for determining the optimal smoothing parameter (p). It is somewhat similar to other model order criteria such as BIC and AIC. Like those, GCV is asymptotically equivalent to "leave one out" cross validation, but requires less computation.

The original GCVSPL code can be run in various modes. The simplest is to directly specify p. Other modes are GCV and "known variance". Oddly, OpenSim only seems to provide the "known variance" mode. So although the tool is named "GCVSpline", the GCV criterion is not supported. Someone please correct me if I am wrong. It would be good if GCVSpline would also support the other modes.

However, after experimenting with the GCV criterion on real data, a long time ago, I concluded that it tends not to smooth enough for using the accelerations. So you have to play around with the smoothing anyway. Also, as John mentioned, and Mohammedreza concluded, it is better to use traditional frequency domain methods where the cutoff frequency is the smoothing parameter. It works just as well and is better understood.

The "known variance" option is attractive for motion capture, because often we know the noise level of the data, e.g. 1.0 mm, where the variance would be 1e-6 m^2. Then the GCVSPL code will find the smoothing parameter p such that the RMS difference between raw and filtered data is exactly the square root of the specified variance.

The smoothing parameter p has a specific meaning. A smoothing spline is the solution to a trajectory optimization problem where the cost function has two terms: data tracking and smoothness. The parameter p is the weighting of the second term. For a cubic smoothing spline, smoothness is the integral of the square of the acceleration. Wikipedia has a good explanation:

https://en.wikipedia.org/wiki/Smoothing_spline.

What is not well known is that a cubic smoothing spline has the same frequency response as a double 2nd order Butterworth filter! I have actually verified this by running artificial sine waves through GCVSPL and Matlab's Butterworth filter. The only difference is in the boundary conditions. A cubic spline assumes zero acceleration at the endpoints. Matlab's filtfilt uses the boundary conditions in [3] (which I have not read). And, as John mentioned, splines can be used with non-equidistant input data.

I have attached the memo that Herman Woltring distributed to GCVSPL users in 1986.

. It provides an equation for the relationship between smoothing parameter p and the cutoff frequency w0 of the equivalent Butterworth filter. My sine wave filtering experiments verified that the equation is correct.

When the GCVSPL code is run in GCV or "known variance" modes, it tells you which optimal value of p it found after the optimization. You can then convert this into an equivalent cutoff frequency using that formula. This would be useful to have available in the OpenSim API.

Ton van den Bogert

[1] Woltring HJ (1986), A FORTRAN package for generalized, cross-validatory spline smoothing and differentiation. Advances in Engineering Software 8(2):104-113 (U.K.).

[2] Woltring HJ (1985) On optimal smoothing and derivative estimation from noisy displacement data in biomechanics. Human Movement Science 4, 229-245.

[3] Gustafsson (1996) Determining the initial states in forward-backward filtering. IEEE Transactions on Signal Processing 44: 988-992.

https://doi.org/10.1109/78.492552.