模型可视化>> gcvspl>> 返回
项目作者： srkuberski

项目描述 ：
GCV spline smoothing
高级语言： C
项目地址: git://github.com/srkuberski/gcvspl.git
创建时间： 2019-11-14T20:59:05Z
项目社区：https://github.com/srkuberski/gcvspl
开源协议：
下载


GCV spline smoothing in Matlab

This repository provides a Matlab interface to Woltring’s classic generalized, cross-validatory (GCV) spline smoothing and differentiation code [1]. Woltring’s original Fortran 77 code [2] was converted to C using the f2c converter [3]. Two Matlab MEX wrappers have been implemented to make the code available to a wider range of users.

Theoretical minimum

  • let $t_i$ ( $i=1,\ldots,n\ge2m$ ) be a set of strictly increasing (not necessarily equidistant) abscissa values with corresponding ordinates $x_i$ and positive weight factors $w_i$
  • a natural spline $s_p(t)$ is a piecewise polynomial function on knot positions $t_i$ which minimizes the criterion function $C_p$ for suitably selected regularization parameter $p\ge0$


Cp=\sum{i=1}^nwi[x_i-s_p(t_i)]^2+p\int{t=-\infty}^{+\infty}|s_p^{(m)}(t)|^2dt

  • through the regularization parameter $p$, a trade-off can be effectuated between the smoothness of the spline (its $m$ th derivative) and the goodness-of-fit to the given data
  • for $p=0$, an exactly interpolating spline is obtained; in the limiting case $p\rightarrow\infty$, the spline becomes an $m$ th order (i.e., an $(m-1)$ th degree) polynomial which fits the data in a weighted least-squares sense
  • Woltring[1] implemented ways to estimate the regularization parameter $p$ by generalized cross-validation or from a given, predicted mean squared error
  • using a basis of B-splines $Bi(t)$, the resulting spline can be written as the linear combination $s_p(t)=\sum{i=1}^nc{p,i}B_i(t)$ with spline coefficients $c{p,i}$
  • the spline polynomials have order $\le2m$ between the knots and order $m$ outside the knot range
  • the polynomials are continuous at the knots up to and including the $(2m-2)$ th derivative with vanishing $m$ th and higher derivatives at the terminal knots $t_1$ and $t_n$
linear cubic quintic heptic
Half order $m$ 1 2 3 4
Polynomial order $2m$ 2 4 6 8
Polynomial degree $2m-1$ 1 3 5 7
Number of continuous derivatives $2m-2$ 0 2 4 6
Minimum number of knots $2m$ 2 4 6 8
Terminal knots w/o full support $2m-2$ 0 2 4 6

Table: Key figures of frequently used spline representations

Matlab interface

For the construction and evaluation of splines, Woltring’s original Fortran code provided two separate functions, gcvspl and splder. Access to these functions is made available here by the use of two Matlab MEX wrappers. In order to use these wrappers in your Matlab installation, you first (and only once) need to compile them on your computer. The compilation requires the presence of a C/C++ compiler software on your computer (e.g., the GNU Compiler Collection gcc). If you do not know if this is the case or not, please refer to the Matlab documentation on how to install, setup, and test a proper MEX building environment before continuing. Once a C compiler is available on your computer, start Matlab, change to the path where the files of this repository are located and run:

  1. mex -v gcvsplmex.c gcvspl.c
  2. mex -v spldermex.c splder.c

Code: Matlab commands to compile the two MEX wrappers

Alternatively, you can use the Makefile provided here. The resulting binaries gcvsplmex.mexa64 and spldermex.mex64 (or similarly named on non-Unix platforms) are Matlab executables which can be used as is. However, this repository also provides (yet another) set of wrappers with additional error handling. The use of these wrappers is recommended. Below you will find short documentation for each of them.

Spline construction with gcvspl

The gcvspl function computes a natural B-spline using the generalized cross-validation and mean-squared prediction error criteria of Craven & Wahba [4]. The model assumes uncorrelated, additive noise and essentially smooth, underlying functions. The independent coordinates may be spaced non-equidistantly.

  1. % compute spline coefficients
  2. %
  3. % c = gcvspl( x, y, m, v, w = ones )
  4. % 
  5. % INPUT
  6. % x : independent variables (numeric row [1, N])
  7. % y : data to be smoothed (numeric matrix [K, N])
  8. % m : spline half order (numeric scalar)
  9. % v : prior given variance [negative: GCV regularization] (numeric scalar)
  10. % w : weight factors (numeric row [1, N])
  11. %
  12. % OUTPUT
  13. % c : spline coefficients (numeric matrix [K, N])
  14. % wk : internal work vector (numeric row [1, 6])
  15. %
  16. % REMARKS
  17. % - array of independent variables x must be sorted
  18. % - spline half orders m = 1, 2, 3, 4 correspond to linear, cubic, quintic, heptic splines (etc.)
  19. % - negative prior given variance v results in generalized, cross-validatory splines
  20. %
  21. % REQUIREMENTS
  22. % - the binary gcvsplmex must be accessible from Matlab's search path

Code: Type help gcvspl to see this in your Matlab installation

Spline evaluation with splder

The splder function computes the values of a B-spline at selected evaluation points. Woltring’s underlying code is an adoption of Lyche et al. [5].

  1. % compute spline values
  2. %
  3. % y = splder( x, c, m, t, n )
  4. %
  5. % INPUT
  6. % x : independent variables (numeric row [1, N])
  7. % c : spline coefficients (numeric matrix [K, N])
  8. % m : spline half order (numeric scalar)
  9. % t : evaluation points (numeric row [1, T])
  10. % n : order of derivative (numeric scalar)
  11. %
  12. % OUTPUT
  13. % y : spline values (numeric matrix [K, T])
  14. %
  15. % REMARKS
  16. % - array of independent variables x must be sorted
  17. % - spline half orders m = 1, 2, 3, 4 correspond to linear, cubic, quintic, heptic splines (etc.)
  18. % - physically, orders of derivate n = 0, 1, 2 correspond to position, velocity, acceleration values (etc.)
  19. %
  20. % REQUIREMENTS
  21. % - the binary spldermex must be accessible from Matlab's search path

Code: Type help splder to see this in your Matlab installation

Zero crossings with splzer

The splzer function computes the location (and sign) of zero crossings of a B-spline.

  1. % compute spline zeros
  2. %
  3. % [z, s] = splzer( x, c, m, n, xmin = min( x ), xmax = max( x ), ofs = 0 )
  4. %
  5. % INPUT
  6. % x : independent variables (numeric row [1, N])
  7. % c : spline coefficients (numeric row [1, N])
  8. % m : spline half order (numeric scalar)
  9. % n : order of derivative (numeric scalar)
  10. % xmin, xmax : search interval (numeric scalar)
  11. % ofs : spline offset (numeric scalar)
  12. %
  13. % OUTPUT
  14. % z : zero locations (numeric row [1, Z])
  15. % s : zero signs (numeric row [1, Z])
  16. %
  17. % REMARKS
  18. % - array of independent variables x must be sorted
  19. % - spline half orders m = 1, 2, 3, 4 correspond to linear, cubic, quintic, heptic splines (etc.)
  20. % - physically, orders of derivate n = 0, 1, 2 correspond to position, velocity, acceleration values (etc.)
  21. % - this function requires the presence of the Curve Fitting Toolbox
  22. %
  23. % REQUIREMENTS
  24. % - Matlab's Curve Fitting Toolbox must be installed on your computer

Code: Type help splzer to see this in your Matlab installation

Matlab example

As an example of usage, consider the following lines of code, which compute the (fifth-order) quintic spline values and their first two derivatives (velocity and acceleration) of a given set of data points (x,y):

  1. c = gcvspl( x, y, 3 ); % compute spline coefficients
  2. y0 = splder( x, c, 3, x, 0 ); % compute spline values
  3. y1 = splder( x, c, 3, x, 1 ); % compute first derivatives (velocity)
  4. y2 = splder( x, c, 3, x, 2 ); % compute second derivatives (acceleration)
  5. z = splzer( x, c, 3, 1 ); % compute zero crossings of the first derivative

Code: A simple example on how to use the splines interface

Note that by specification of more (or less) than x evaluation points in the fourth argument to the function splder, re-sampling of the original data is possible.

License

For archiving and documentation purposes, the original Fortran 77 source file (gcvspl.f) is duplicated here, as well as the intermediary C code files (gcvspl.c, gcvspl.h). Please adhere to their original copyright (unrestricted non-commercial use) [2]. Other software provided with this repository is released to the public domain.

If you use this software in your work, please cite it using the following metadata:

Kuberski, Stephan R. (2023). GCV spline smoothing in Matlab [Software]. url: github.com/srkuberski/gcvspl

References

[1] Woltring, Herman J. (1986). A Fortran package for generalized, cross-validatory spline smoothing and differentiation. Advances in Engineering Software 8(2). doi: 10.1016/0141-1195(86)90098-790098-7).

[2] International Society of Biomechanics. Signal processing software. url: isbweb.org/software/sigproc.html. (retrieved in November 2019).

[3] AT&T Bell Laboratories. A Fortran to C converter. url: www.netlib.org/f2c. (retrieved in November 2019).

[4] Craven, Peter and Grace Wahba (1979). Smoothing noisy data with spline functions. Numerische Mathematik 31(4). doi: 10.1007/BF01404567

[5] Lyche, Tom, Larry L. Schumaker, and Kamy Sephehrnoori (1983). Fortran subroutines for computing smoothing and interpolating natural splines. Advances in Engineering Software 5(1). doi:10.1016/0141-1195(83)90073-690073-6)