C*********************************************************************** C C GCVTST.FOR, 1986-02-11 (comm. 1988-11-15; logical corr. 1991-09-03) C C Author: H.J. Woltring C C Organizations: University of Nijmegen, and C Philips Medical Systems, Eindhoven C (The Netherlands) C C----------------------------------------------------------------------- C C Testprogramme for generalized cross-validatory spline smoothing C with subroutine GCVSPL and function SPLDER using the data of C C.L. Vaughan, Smoothing and differentiation of displacement- C time data: an application of splines and digital filtering. C International Journal of Bio-Medical Computing 13(1982)375-382. C C These data describe a constant acceleration paradigm (i.e., a C falling golf ball). Therefore, the true third and higher deri- C vatives are zero, so constraining is appropriate for the third C derivative only, i.e., a quintic spline. For lower orders, an C erroneous constraint is imposed on acceleration, velocity, or C displacement; for higher orders, acceleration may become some C arbitrary polynomial in such a way that the displacement data C exhibit a reasonable goodness-of-fit. For other kinds of sig- C nal, constraining one derivative higher than the highest one C sought is the appropriate choice. C C The variable FRE indicates the approximate cut-off frequency of C the spline's equivalent Butterworth-filter (in Hz). For equi- C distantly sampled, uniformly weighted data, this filter has a C transfer function [ 1 + (F/FRE)^2M ]^-1, where F is the fre- C quency, and M the spline's order. Effectively, the spline be- C haves as two cascaded, M-th order Butterworth filters without C phase distortion. C C The programme types out statistics on the estimation procedure C and on the estimated second derivatives. If the DEBUG-lines are C compiled, also the raw data and the estimated spline values, C first and second derivatives at the knot positions are typed. C C The only subprogrammes to be conventionally called by a user C are subroutine GCVSPL for calculating the spline parameters, C and function SPLDER for calculating the spline function and its C derivatives within the knot range. See the comments in the C headers of these subprogrammes for further details. C C The TYPE and ACCEPT statements are typical for VAX/VMS; if so C needed, they can replace the WRITE and READ statements in the C subsequent lines. These interactive statements are present in C the main programme only; the subroutine package should be fully C portable under FORTRAN-77. C C Reference: C C H.J. Woltring (1986), A Fortran package for generalized, cross- C validatory spline smoothing and differentiation. Advances in C Engineering software 8(1986)2, pp. 104-113. C C********************************************************************** C PROGRAM GCVTST !REAL*8 C IMPLICIT REAL*8 (A-H,O-Z) PARAMETER ( K=1, NN=50, MM=10, MM2=MM*2, NWK=NN+6*(NN*MM+1) ) PARAMETER ( PI=3.1415926536, SCALE=0.5/PI, lout=6, lin=5 ) DIMENSION X(NN), Y(NN), WX(NN), C(NN), WK(NWK), Q(0:MM), V(MM2) C DATA Y/1.770D0,1.757D0,1.748D0,1.740D0,1.726D0, !Noisy samples 1 1.715D0,1.698D0,1.683D0,1.667D0,1.651D0, !of the vertical 2 1.632D0,1.612D0,1.593D0,1.572D0,1.551D0, !coordinates (in 3 1.530D0,1.507D0,1.483D0,1.445D0,1.428D0, !metres) of a 4 1.401D0,1.371D0,1.343D0,1.311D0,1.279D0, !freely falling 5 1.245D0,1.212D0,1.175D0,1.143D0,1.105D0, !golf ball. 6 1.063D0,1.029D0,0.991D0,0.953D0,0.910D0, 7 0.869D0,0.823D0,0.779D0,0.732D0,0.691D0, 8 0.644D0,0.595D0,0.548D0,0.501D0,0.447D0, 9 0.395D0,0.350D0,0.294D0,0.243D0,0.185D0/ DATA WX/50*1D0/, WY/1D0/, AT/9.85D-3/ !Weights, sampling interval C C*** Create time array (knot array) C DO 10 IX=1,NN X(IX) = AT * (IX - 1) 10 CONTINUE C C*** Get parameters (see comments in subroutine GCVSPL) or exit C C M: spline order (0=exit, C 1=linear, 2=cubic, 3=quintic, C 4=heptic, 5=nonic, ...) C C |MODE|: 1=non-iterative, fixed p-value C 2=GCV C 3=variance C 4=degrees-of-freedom C MODE < 0: retry for same M, |MODE|, N C C VAL: value if appropriate for selected MODE C C N: number of points to be used ( N <= 50 ), default 50. C C*** C C 20 TYPE 700 20 write(lout,700) 700 FORMAT(/'$M, MODE, VAL , N = ') C ACCEPT 710, M, MODE, VAL, N read(lin,710) m, mode, val, n 710 FORMAT(2I10,E15.0,I10) IF ((N.LT.2*M).OR.(N.GT.NN)) N = NN IF ((MODE.EQ.0).OR.(IABS(MODE).GT.5).OR. 1 ( M.LE.0).OR.( M.GT.MM)) STOP 'GCVTST: ready' !exit C C*** Assess spline coefficients and type resulting statistics C CALL GCVSPL ( X, Y, NN, WX, WY, M, N, K, MODE, VAL, C, NN, 1 WK, IER) IF (IER.NE.0) THEN C TYPE 720, IER write(lout,720) ier 720 FORMAT(' GCVTST: error #',I3) GO TO 20 !next trial ELSE VAR = WK(6) IF (WK(4).EQ.0D0) THEN FRE = 5D-1 / AT ELSE FRE = SCALE * (WK(4)*AT)**(-0.5/M) ENDIF C TYPE 730, VAR, (WK(I),I=1,4), FRE write(lout,730) var, (wk(i),i=1,4), fre 730 FORMAT(' var =',1PD15.6,', GCV =',D15.6,', msr =',D15.6/ 1 ' df =',0PF8.3,', p =',1PD15.6, 2 ', fre =',1PD15.6) ENDIF C C*** Reconstruct data, type i, x(i), y(i), s(i), s'(i), s''(i) [D] C*** Assess and type acceleration mean and standard deviation C C TYPE 740 CD write(lout,740) CD740 FORMAT('0 i x(i) y(i) ', CD 1 ' s0(i) s1(i) s2(i)'/) IDM = MIN0(2,M) DACCAV = 0D0 DACCSD = 0D0 Q(2) = 0D0 DO 40 I=1,N J = I DO 30 IDER=0,IDM Q(IDER) = SPLDER ( IDER, M, N, X(I), X, C, J, V ) 30 CONTINUE C TYPE 750, I, X(I), Y(I), (Q(IDER), IDER=0,IDM) CD write(lout,750) i, x(i), y(i), (q(ider), ider=0,idm) CD750 FORMAT(I4,':',5F13.6) DACCAV = DACCAV + Q(2) DACCSD = DACCSD + Q(2)*Q(2) 40 CONTINUE ACCAV = DACCAV / N ACCSD = DSQRT((DACCSD - ACCAV*DACCAV)/(N-1)) C TYPE 760, ACCAV, ACCSD write(lout,760) accav, accsd 760 FORMAT('0acceleration mean and sd:',2F12.7,' m/s**2') C C*** Next run C GO TO 20 END