1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
|
# Copyright(c) 1986 Association of Universities for Research in Astronomy Inc.
include "bspln.h"
.help spllsq 2 "math library"
.ih ___________________________________________________________________________
NAME
spllsq -- general smoothing b-spline with uniform knots.
.ih
USAGE
spllsq (x, y, w, npts, bspln, q, work, k, n, wflg, ier)
.ih
PARAMETERS
See the listing of SPLSQV if a more detailed discussion of the parameters
than that given here is desired.
.ls x,y,w
(real[npts]). Abscissae, ordinates, and weights of the data points.
X[1] and X[NPTS] determine the range of the spline (the interval over which
it can be evaluated). Dummy data points with zero weight can be supplied
to fix the range if desired.
.le
.ls npts
The number of data points.
.le
.ls bspln
(real[2*n+30]). Output array of N b-spline coefficients, N+K knots,
and the spline descriptor structure.
.le
.ls q
(real[n,k]). Work array.
.le
.ls work
(real[n]). Work array.
.le
.ls k
The order of the desired spline (1-20) (Cubic == 4).
.le
.ls n
N is the number of b-spline coefficients to be output. The relationship
between N and the number of polynomial pieces in the spline (NPP)
is given by NPP = N-K-1. Note that NPP is the important parameter.
N is used in the parameter list rather than NPP because it is used in
dimensioning the arrays.
.le
.ls wflg
Weight generation flag. Options:
.nf
0 pass W on to SPLSQV as is
1 calculate default weights
2 same as 1, but set W[i] only
if W[i] already NONZERO
.fi
If WFLG==2, data points may be rejected from the fit by setting their
corresponding weight to zero, before calling SPLLSQ (good points must be
assigned nonzero ws before calling SPLLSQ). Use WFLG==1 if no points
are to be rejected, to avoid having to initialize the W array on every
call.
.le
.ls ier
Error return. Zero if no error. Nonzero indicates invalid combination
of N, K, and NPTS.
.le
.ih
DESCRIPTION
A general algorithm for smoothing with uniform B-splines of arbitrary order.
Calls SPLSQV, which is adapted from L2APPR, as described in "A Practical
Guide To Splines", by C. DeBoor.
SPLLSQ is identical to SPLSQV, except that (1): the array T[n+k], giving
the x-positions of the knots of the b-spline, is generated internally,
rather than by the calling program, and (2): the weights associated with
the data points may be calculated automatically. The knots are generated
with a uniform spacing. Note that the x-values of the data points do not
have to be uniformly spaced, but they must be monotonically increasing, and
(1) both must span the same range, and (2) if N spline coefficients are
desired, there must be at least N+K data points.
.ih
BUGS
A routine FSPLSQ needs to be written to make use of the Cholesky factorization
returned by SPLSQV in W1, for more efficient fitting of subsequent data sets
that differ only in the y-values of the data points.
.ih
SEE ALSO
seval(2)
.endhelp _______________________________________________________________________
procedure spllsq (x, y, w, npts, bspln, q, work, k, n, wflg, ier)
real x[npts], y[npts], w[npts]
real q[n,k], work[n]
real bspln[ARB]
int wflg, npts, n, k, ier
int i, npp, km1, knot
real dx
begin
ier = 0 #successful return value
npp = n - k + 1 #number of polynomial pieces
if (npp <= 0) {
ier = 1
return
}
# Set up spline descriptor in the array BCOEF, for later evaluation
# of the spline by BSPLN (see "bcoef.h" for definitions of the fields).
NCOEF = n #number of b-spline coeff
ORDER = k #order of the spline
XMIN = x[1] #x at left endpoint
XMAX = x[npts] #x at right endpoint
KINDEX = KNOT1 - 1 + k #for SEVAL
# Set values of knots. First K knots take on the value of the first
# breakpoint, next npp-1 knots have spacing DX, and the last K knots
# take on the value of the last breakpoint.
dx = (XMAX - XMIN ) / npp #span(x) = span(t)
km1 = k - 1
knot = KNOT1 - 1
do i = 1, km1
bspln[knot+i] = XMIN
knot = knot + km1
do i = 1, npp
bspln[knot+i] = XMIN + (i-1) * dx
knot = knot + npp
do i = 1, k
bspln[knot+i] = XMAX
# Calculate default weights. The default weight of a datapoint is
# proportional to the spacing.
if (wflg > 0) {
if (w[1] > 0 || wflg == 1)
w[1] = x[2] - x[1]
do i = 2, npts - 1 {
if (w[i] > 0 || wflg == 1)
w[i] = (x[i+1] - x[i-1]) / 2.
}
if (w[npts] > 0 || wflg == 1)
w[npts] = x[npts] - x[npts-1]
}
# Call SPLSQV to do the actual spline fit. The N b-spline coefficients
# of order K, N+K knots, and the spline descriptor are returned in
# BSPLN, for evaluation of the spline via SEVAL.
call splsqv (x, y, w, npts, bspln[nint(KNOT1)], n, k, q,
work, bspln[COEF1], ier)
end
|
