LAMONT DOHERTY GEOLOGICAL OBSERVATORY
OF COLUMBIA UNIVERSITY
TREE-RING LABORATORY Palisades, New York 10964
*
Technical Report No. CU-1-81/TR1
This material is based upon work supported by the National Science Foundation under Grant No. ATM79-19223, Gordon C. Jacoby, Principal Investigator.
*
THE CUBIC SMOOTHING SPLINE AS A DIGITAL FILTER
by
Kenneth Peters and Edward R. Cook October 1981
*
Reproduction in whole or in part is permitted for any purpose of the United States Government. Others must secure the author's permission for use of this material.
★
Distribution of this document is unlimited.
LAMONT- DOHERTY GEOLOGICAL OBSERVATORY of Columbia University TREE-RING LABORATORY Palisades, New York 10964
Technical Report No. CU-1-81/TR1
*
This material is based upon work supported by the National Science Foundation under Grant No. ATM79-19223, Gordon C. Jacoby, Principal Investigator.
*
THE CUBIC SMOOTHING SPLINE AS A DIGITAL FILTER
by
Kenneth Peters and Edward R. Cook October 1981
*
Reproduction in whole or in part is permitted for any purpose of the United States Government. Others must secure the author's permission for use of this material.
*
Distribution of this document is unlimited.
■k
Any opinions, findings, and conclusions or recommendations expressed in this publication are those of the authors and do not necessarily reflect the views of the National Science Foundation.
Digitized by the Internet Archive in 2020 with funding from Columbia University Libraries
https://archive.org/details/cubicsmoothingspOOpete
THE CUBIC SMOOTHING SPLINE AS A DIGITAL FILTER
Kenneth Peters and Edward R. Cook
ABSTRACT
Reinsch (1967) introduced the cubic smoothing spline as the solu¬ tion to a certain variational problem in the theory of data smoothing.
This smoothing operator has good transient response, a minimal curva¬ ture property, and compensation for end effects so that all the data is used. In addition, the data can have arbitrary spacing and weighting.
In this paper we represent Reinsch1 s equations in terms of convolu¬ tions for the special case of equally-spaced and equally-weighted data.
We then use a Fourier transform to obtain a frequency domain represen¬ tation, solve the characteristic polynomial of the system, and apply an inverse transform to the frequency response function to obtain the impulse response function. The Lagrangian multiplier of the variational problem parameterizes the family of digital filters derived in this man¬ ner. This representation which is exact for finite data sets extends the range of application of the smoothing spline and, incidentally, simplifies the algebra and reduces the amount of computer time and storage required, by allowing one to select the degree of smoothing on the basis of frequency domain considerations. One can do time domain filtering without end effects in a computationally efficient manner.
Smoothings of four tree-ring width sequences from forest interior sites are presented to illustrate the particular suitability of these filters for data exhibiting both non-stationary and episodic behavior.
1
THE CUBIC SMOOTHING SPLINE AS A DIGITAL FILTER
I INTRODUCTION
Reinsch [7] discusses the problem of smoothing noisy data when the underlying function is not known. We begin with a brief statement of his main results in his own notation and terminology. He suggests that the fitting function, g(x), minimize the total square curvature,
2
dx
(la)
xonder the constraint that a weighted sum of squares of the residuals be bounded above ,
n
"g(Xi)-yil 2 < S;
6Yi
(lb)
where the y^ represent the noisy data, the 6y^ are a set of weights and S is a scaling parameter. Equivalently, introducing a Lagrangian multiplier, p, and an auxiliary variable, z, the functional
is to be minimized.
The solution is a cubic spline, i.e., a concatenation of cubic segments :
g^(x) = ai + (x-x^) + c^ (x-x^) 2 + d^ (x-x^) 3, x^ ^ x < x^+1, (2)
such that g(x) , dg(x)/dx, and d2g(x)/dx- are continuous at x^ for i=l, n-1.
2
If the Lagrangian multiplier, p, is known, the values can be computed from ^
(QD Q + pT ) c = pQy, (3)
where the circumflex denotes transposition, and the a^ from
Then
a = y -
p-1D2Qc.
(4)
di = (ci+1 - ci) / 3hi, i = 0, n - 1 (5)
where, for the sake of uniformity of notation, Reinsch chooses to set Cq and c^ equal to 0 ; and
bi = E(ai+1 ~ a±) /hi3 - cihi - dihi , i = 0, n - 1. The notation used is Reinsch' s:
hi = *
i+1 " xi;
c (c]_, • • • 'cn-l^ '
(6)
(7a) (7b) (7c) ( 7d) (7e) (7f)
y = (y0# yi . yn.r yn) *
a ^ a0 ' al' • • • '^-1' an^ '
D = diag (6yQ, . . . ,6yn) ;
T = a positive-definite tridiagonal matrix of order n-1,
where t^ = 2 (hi_1 + hi)/3, and ti>i+1 = h±/3 = ti+1/i;
Q = a tridiagonal matrix with n+1 rows and n-1 columns (7g)
where qi_1^i = l/hi_1, qi;j_ = - (l/hj^ + l/hi) , and
gi+l,i " 1/bi*
For Reinsch1 s intended application to data smoothing it is natural to specify the weights 5y^, usually as the square root of the error variance at x=x^ if it is known or can be estimated, and S, as the value
of a chi-square variable with n degrees of freedom. He then proceeds to
3
and may be of interest to the reader. Reinsch t 8] recasts the problem of reference £71 in much more general terms and shows that similar results hold. Demmler and Reinsch £41 achieve by advanced methods a very general result of which formula (16) of this paper is a special case. Wahba £9] discusses the problem of optimizing the degree of smoothing given certain assumptions regarding the random component of the data and gives explicit formulas which can be used in conjunction with the heuristic considerations which motivated this paper. Wahba and Wold £10] and Craven and Wahba [3] apply a method of cross-validation to select the smoothing parameter and in £11] Wahba and Wold discuss the use of periodic splines in spectral density estimation which is related to a suggestion for natural splines in this paper in Section V. Wold £12] is concerned with the suitability of splines of various kinds from a data analytic point of view. In this paper we are concerned with a discrete problem — the relationship between the vectors y and a. In the limit of large p there are no digital filtering effects. Horowitz C6] considers the continuous problem of the filtering effects of common interpolating splines. His analysis may be of interest if the dis¬ crete data arise from accumulation or from sampling continuous records.
The references given here should provide an adequate introduction to smoothing splines.
The rest of the paper consists of four sections and two appendices. Section II contains the derivation of the convolutional representation and of the frequency response function. Section III contains the solution of the characteristic equation of the system. Section IV contains the deri¬ vation of the impulse response function. Section V includes a discussion of how to interpret and apply the representation derived in sections II,
III and IV, and an example which illustrates the time domain behavior of
5
the cubic smoothing spline. Appendix A is a discussion of this behavior in
the limit of small curvature, and Appendix B contains an analysis of the coupled ordinary differential equations of which the cubic smoothing spline is an interpolated finite difference solution.
The reader interested only in practical applications should read sections I and V and consider the main practical results which are in equations (16) and (33) and figures 1 and 5. The critical reader will want to include the beginning of section II. The appendices and the parts of the text relating to the correction terms (indicated by primes) and the small curvature limit, can be ignored, but clues to understanding the Fourier representation are contained therein. Finally, those filters generated by the real negative roots of the characteristic polynomial are of little practical interest since they result in very little smoothing.
II THE CONVOLUTION REPRESENTATION AND THE FREQUENCY RESPONSE FUNCTIONS
Consider the special case of equations (3) and (4) for which h^ = 1.0 = Sy^ for all i. In this case the matrices Q and T and R =
QQ + pT will be Toeplitz, i.e. , all the terms on the same diagonal will be equal, and symmetric. The matrix- vector product Rc will be equivalent to convolution with a symmetric (phase-free) filter except for the first and the last two points. In fact these restrictions are necessary because convolution is defined as a time-invariant linear process, and it is impractical to do Fourier transforms on unequally spaced data.
The particular numerical values of the spacing, in this case unity,
-3
serves to scale or normalize the parameter p which has the units of x With these restrictions then the desired representation can be achieved as follows.
6
First, extend a typical row of the matrix R by appending an infinite number of zeroes on each end, and do the same for Q and Q. This exten¬ sion is necessary because the limits on the Fourier integral are infinite. The particular form of the extension (zeroes) is dictated by the band form of the corresponding matrices. Thus define
r = {...,0,1, (-4+p/3) , (6+4p/3) , (-4+p/3) ,1,0, . . . } (8)
and
q = {... ,0,1, -2, 1,0,.. .} (9)
in accordance with (7f) and (7g) . The subscript p makes the para¬ meterization explicit. The vectors c, y and a can be extended in a similar manner, i.e.,
c = {.. . ,0,cQ,c1,... ,cn_1,cn,0, .. .}, c0=cn=0, (7b')
y = {...,0,y0,y1,...,yn_1,yn,0,...}, (7c1)
and
a= { . . . ,0,aQ,a1, . . . ,an_1,an,0, . . . }; (/a')
but this particular extension for c, y and a is not dictated by
Reinsch's equations (except for the end conditions Cg=cn=0) which
entail no assumptions regarding the behavior of y outside the data
interval. We use it only for definiteness and convenience. More
general extensions will be discussed in section V. The main point
here is to represent Reinsch's computations as exactly as possible.
Next notice that while (3) equates two vectors of length n-1,
r * c, where * denotes convolution, and q * y will in general have n + 3 J?
non-zero elements. To represent (3) accurately the extra terms gener¬ ated by the convolutions must be subtracted from both sides before they are set equal. Thus, (3) becomes
rp*c - = p (q * y-y') ( 3 ’ )
7
with
c ' = { . . . ,0,0^ , i (-4+p/3) c^+c^ ,0 , . . . ,0, Ccn_2+ (-4+p/3) cn_12 , cn_]_,0, . . . } (10)
and
y' = { • • • /0,y0, (-2y0+y1) ,0, . . . ,0, (-2yn_1+yn) ,yn,0, . . . }. (11)
Note that the prime will not be used to denote differentiation in this paper.
In the case of equation (4) the extra terms are retained. It becomes
y - — q P
c, p > 0,
(4')
Finally, to achieve symmetry in the transforms, and for representa¬ tional purposes, the discrete vectors will be represented as weighted sums of delta functions. Letting a denote any of them define
a(t) 7 a. 6(t-i). (12)
*• — l
i=-°°
The explicit dependence on t will identify continuous time domain functions.
Fourier transforming (3') and (4') modified according to (12)
|
gives |
r (f)c(f) - Cp ( f ) = p Cq(f)y(f) - y'(f)] |
(13) |
|
|
and |
a (f ) = y(f) - |q(f)c(f) , p > 0. ir |
(14) |
|
|
Solving |
(13) |
for c(f), substituting in (14) and rearranging gives |
|
|
a (f) |
= u (f)y(f) + u ' (f) Cy' (f) - ic • (f)] , p > 0, p P P P |
(15) |
|
|
where |
u (f) |
= (p/6) cos2TTf+ (p/3) |
(16) |
|
P |
cos^2'n'f+ ( (p/6) -2) cos2nf+ ( (p/3) +1) |
||
|
and |
n , _ 2 (cos2TTf-l) |
(17) |
|
|
cos^2iTf+ ( (p/6) -2) cos27rf+ ( (p/3) +1 |
|||
|
which are the |
frequency response functions. |
||
|
The |
functions u (f) and u '(f) are shown in figures 1 'and 2 over P P |
the |
|
|
useful |
range |
of p. |
8
Ill THE ROOTS OF THE CHARACTERISTIC POLYNOMIAL
The characteristic polynomial of the system is the (common) denomi¬ nator of u (f) or u ' (f) expressed in terms of z = exp(27rif) . The roots I? P
of this function will be needed in section IV. Setting the denominator of Up(f) equal to zero and solving for cos2irf gives
and
Zp = 1- (p/12) + i(p(72-p))1/2/12
z" = 1- (p/12) - i (p (72-p) ) 1/2/12
(18)
(19)
The principle root is used here and always. These roots lie on the
circle of radius 3 centered on the real axis at -2 for 0<p$72.
For p£72 Im z =0. Also, to avoid ambiguity, lim z as p-*30 is defined P P
to be -2.
If z=exp (27rif ) is a root, so is 1/z since cos 2nf = (exp(2irif) + exp (-27rif ) ) /2 . Using this relationship and its converse, the four roots are, for 0<p<72,
and
++ + • +2» 1/2
z = z + i ( 1- z )
P P P
H — + . , +2 . 1/ 2
z = z - i(l-z )
p p p
-+ - -2 1/2 Zp = zp + i(l-zp )
-2 1/2
z = z - i(lz )
P P P
(20)
(21)
(22)
(23)
For p^72 there are four distinct negative real roots. Two of
these, z++ and z +, remain within the unit circle,
P P
and
z + P
(24)
(25)
9
while the reciprocals, z^+ and z ^ , respectively, remain outside. The
lim z as p-*°° is defined to be 0, while lim z as p-*» is -2 + /3.
P P
The root z^"*" is plotted in figure 3. The other roots are simply related to z^+. Thus z^ is the reciprocal, z^~ is the complex con¬ jugate, and z is the reciprocal of the complex conjugate.
P
IV THE IMPULSE RESPONSE FUNCTIONS
Using the results of section III, the impulse response functions
can be obtained as the inverse Fourier transforms of the frequency
response functions u (f) and u ' (f) . Decomposing the right side of (16)
P P
into partial fractions gives
u (f) =
Ar
Br
cos2,n'f - z2
cos2iTf - zT
(26)
where
A =
(p/6) Zp + (p/3)
z+ - Z~
and Br
(p/6) Zp + (p/3)
z~ — z+ p p
(27)
Then
u (t) = A ! exp(2TTift) df + B j exP (2irif t)
COs2fTf - z"
P t cos2irf - zx
d p
df
(28)
P
Tri
z^dz
z2-2z+ + 1
P
+ !p
ttT
z^dz
oo
E
5 (t-i)
^ z -2z z + 1
<* p
C *
(29)
the unit circle, C, being traversed in the positive sense for t>0. To see that (29) follows from (28) consider (12) and the fact that the series of delta functions used to repeat the principal part of u^(f) by convolution is the transform of the series in t.
10
The integrals in (29) can be evaluated by the method of residues Cl] .
The roots of the denominator of the first integrand are z++ and z+~ ,
P P
and for the second z and z , which are respectively complex conjugates
P P
of the first (reciprocal) pair. Then, for all p>0, the impulse response function
u (t)
P
(2K z++|t|+ 2 P P
V
-+ t
CO
) .XI <5 (t-i) =
l=-oo
V (t)
P
.£
l=_oc
6 (t-i) ,
where
K =
++
A / (z P P
++
Vzp ) ,
and
L„ =
B /(z P P
-+
Vzp ) .
(30)
(31)
(32)
Using the same conjugacy relationships between the roots, for 0<p<72 ,
Vp(t) = 4Re{Kpz++ I t I } = 4|Kp| |Zp+| I tlcos (Arg Kp + |t|Arg Zp+) (33)
In the limit as p-K), vp(t) -> 0 uniformly everywhere.
For p>72 there are four distinct negative roots occurring in reciprocal pairs. The impulse response function becomes the sum of two complex numbers with different amplitudes (decaying exponentially in |t|) rotating in opposite senses around the t axis, piercing the imaginary-t plane at half integral values of t and the real-t plane at integral values of t. In the limit as p-*00, v (t) 1 for t=0 and -*■ 0 elsewhere. That is, u (t) is a
ir ir
spike at the origin.
The function v^ (t) must be evaluated by a limit process. The limits from the right and the left are the same and thus Vp(t) is continuous at p=72 for all values of t. The function v (0) in particular is plotted in figure 4 .
11
The functions v^(t) for log p = -3, -4, -5
The analysis for u ' proceeds as for u
P p
and B ' are different.
P
are plotted in figure 5 . except that the constants
A ' = (2z+ - 2) / (z+ - z") P P P P
and
(34)
-3
V = (2zp ' 2,/(zp " ZP (35)
v ' (0) is plotted in figure 4, and the functions v ' (t) for log p = P P
-4,-5 are plotted in figure 6.
V DISCUSSION AND AN EXAMPLE
Transforming (15) back into the time domain and reverting to the discrete vector symbolism gives
a = u*y-u'* (y ' - — c' ) , p > 0, (36)
1 p 7 p p ^
and the operation of the smoothing spline is seen to be representable in the time domain as the sum of two convolutions. The first term accurately represents the operation when end effects can be ignored; for instance far enough away from the ends of the data so that the filter cannot feel them, or in the cases in which a reasonable extra¬ polation of the data can be made. The second term is peculiar to the smoothing spline and accounts for its apparent lack of end effects.
It is necessary to consider the sense in which the smoothing spline has no end effects.
This phrase is appropriate if one accepts the global variational definition of the smoothing spline as implying that its quasi-local behavior is the same at the ends of the data set as near the center.
12
Indeed the Euler-Lagrange equations must be satisfied at all interior
points. In this sense the second term in (36) "fixes" the computations
so that u alone is an adequate description of its operation, in which P
case one need only consult figure 1 in selecting an appropriate p value. We have been using the spline in this sense. It is more appropriate in this paper, however, to consider the end behavior in terms of the con¬ volution representation.
If a different choice had been made in extending the vectors c, y and a, the contribution of the second term in (36) would have been dif¬ ferent. In fact there are (an infinite number of) extensions such that this second term would contribute nothing. They can be computed by con¬ volving Up with y and requiring that the resulting values be equal to a^ for 0<i<n. If one accepts these extensions as "reasonable", then again u^ by itself adequately describes the filter. This implies that the smoothing spline can be used as a predictor or to achieve higher reso¬ lution in spectral estimates. If the physical model entails that the data vary according to some principle of minimum curvature at least at lower frequencies such an idea is certainly tenable, but as in other kinds of prediction and spectral enhancement (maximum entropy, for instance) due consideration must be given to the underlying model. Prediction is reliable only if the model is appropriate.
Four smoothings of tree-ring width sequences are shown in figure 7. The same figure appears in [2] . The dashed lines show orthogonal poly¬ nomial fits generated by the INDXA program C53 the orders of which are determined by stopping when an additional order fails to reduce the residual variance by 5% or more. The log p value was -4 in each case.
13
Tree-ring width sequences from forest interior sites tend to exhibit behavior similar to that seen in records of complex physical processes such as level changes in streams and reservoirs and climatic variations. These data exhibit sudden changes in level (episodic behavior) over a wide range of time scales. The spectra tend to be lumpy and have concen¬ trations of power near zero frequency. The particular suitability of the smoothing spline as a low pass filter for certain tree-ring records, and perhaps for other similar data, is accounted for by the shape of the impulse response functions. The spikey exponential envelope with the |t| dependence ensures that the filter will be responsive to sudden changes in level. On the other hand, the fact that it is a low-pass filter with a minimum curvature property ensures that power near zero frequency will be preserved. The frequency response function is smooth at the origin. Thus trends and slow wandering in the data will be preserved. Selecting a value of S (or p) affects a compromise between fitting and smoothing the data. We see how this compromise manifests itself in the digital filters .
Programming of equations (31) and (4') can be done in various ways.
We have not tried to optimize the computations beyond using an equation solver specialized for symmetric band matrices. The user concerned about storage might try a Toeplitz-specif ic algorithm. The only caution regards the possible necessity of using double precision accumulation variables.
We found that in single precision (32 binary bits) the computations were
-4
unstable for values of p smaller than 10 roughly.
14
-
APPENDIX A
It is not obvious that the weight function shown in figure 6 will, in the limit of small p, give a series of values which lie on a straight line for O^i^n and be zero outside this interval when convolved with anything at all. The purpose of this appendix is to examine this limit¬ ing case.
For 0<p<72, and thus in the limit as p -*■ 0 ,
II 1 1 1
v ' (t) = 4Re{Kp ' z^+ ' 1 }=4 | Kp ' | | z++ | cos (Arg Kp+|t|Arg z^+) . (Al)
Letting r, 0' and 0 denote small positive numbers, for small p Al can be written
v ' (t) = 4 Ik ' | (l-r) 1 1 1 cos ( (3tt/4) -0' - 1 1 1 0 ) . (A2)
p P
Expanding in r and 0 ' and 0 ,
v 1 (t) = k [1 - 1 1 1 r + o(r2)] El-0' - 1 1 1 0 - o(02)], (A3)
P
where
k = -4|K ' |/V2.
To first order therefore
v ' (t) = k(l-0‘- s 1 1 1 ) , for s < < 1 1 1 < < 1/s, (A4)
ir
where s = r + 0.
1
Denoting the four non-zero values of (y'-vc ') by w ,
ir P “ 1
wn+l' at least the following four requirements must be met:
k(l-0,)w_1 + k (1-0 ' -s) wQ = 0,
k(l-@'-s)w_1 + k (1-0 ' ) wQ = aQ,
k(l-0')wn + k(l-0'-s)wn+1 = an,
and k(l-0'-s)w + k(l-0')w = 0.
n n+1
w
O'
w
n
(A5)
(A6)
(A7)
(A8)
15
Since in the limit of small p,|K 1 | and thus |k| become very large,
ir
k can be divided out of A5 - A8. Solving the resulting equations for
w and w (A5 and A6) and w and w (A7 and A8) yields to first order -1 0 n n+1
in s ,
|
w_x = -aQ (1-0 ' -s) /2s ( 1-0 ' ) , |
(A9) |
|
|
wQ = a0/2s. |
(A10) |
|
|
w - a /2s, n nx ' |
(All) |
|
|
and |
wn+l = -an(l-0,-s)/2s(l-0') . |
(A12) |
Since s > 0 for all finite p these solutions are defined for the limit. Also, by the symmetric linear (to first order) dependence of v ' (t)
IT
on t, it follows from A5 that a^ = 0 for all -1/s <<i< -1; from A8, that a^ = 0 for all n+l<i<<l/s; and from A6 and A7 that a^ varies linearly with i (to first order) for O^i^n. In the limit as p-*0 the first order expressions become exact and l/s-*-°°.
16
APPENDIX B
The first terms in the central finite difference expansions of the second and fourth derivative operators are obvious in the matrices Q and R. Corresponding to equations (3') and (4') there are the coupled differential equations
CD4 + (p/3) (D2 + 6)] C(t) = pD2Y(t) (Bl)
and
A ( t) = y - — D2C ( t) , p>0 . (B2)
P
Conversely, a finite difference solution of (Bl) and (B2) will satisfy equations (3') and (4'). More generally, the theory of Reinsch's smoothing splines is related to the theory of the interpolated finite difference solutions of certain ordinary differential equations with constant coefficients. The remainder of this appendix is a brief analysis of the equations (Bl) and (B2) , analogous to that done for the spline equations, which will facilitate some comparisons the reader may wish to make between the continuous and discrete systems.
Fourier transforming equations (Bl) and (B2) and solving for A(f) in terms of Y(f) , gives
p(3-27T2f2)
A (f ) = - — - — — Y(f) =U (f)Y(f). (B3)
24irf +p(3-2tt f2)
Note that this same transform would arise from the single equation
(D4 + (p/3) (D2 + 6))A(t) = (p/3) (P2 + 6) Y ( t) .
(B4)
This won't work in the discrete case because of the correction terms. The roots of the denominator for 0$p-$72 are
fp+ = (P + i (p (72-p) ) 1/2 ) /24tt2
cin d i ^ p
f2_ = (p - i (p (72-p) ) )/24tt
ir
(B5)
(B6)
17
and therefore
and
+<£p+)1/2-
-if2*) 1/2
P
-(f2_)1/2.
P
(B7)
(B8)
(B9)
(BIO)
Corresponding to the constants A and B
P P
there are
2 2+ 4 9+ o_
a = p (3-2tt f ) /24tt ( r + - f ), P P P P
and
2 2- 4 2- 2+
3 = p(3-2rr f ) /24ir (f - f ) ,
P P ir P
in the inverse transform
(Bll)
(B12)
while
k
P
Vf
++ P '
(B13)
and
£P - vf;+
(B14)
The contour of integration consists of the real axis and the upper half of the point at infinity, traversed counterclockwise for t>0.
For 0<p<72
U (t) = -2iTexp (-2iTlm{ f^+} 1 1 1 ) Im{k^exp (27riRe{fp+} 1 1 1 ) } . (B15)
For p^72, each contour (t<0) includes two of the real roots and
U (t) = -ir(k sin 2tt f 1 1 1 + 1 sin 2ir f_+ 1 1 1 ) . (B16)
P P P P P
There is no decay in the continuous case for p^72.
Evaluation of the limits as p approaches 72 from the left and right
will show that for any value of t, U (t) is continuous at p=72.
P
18
ACKNOWLEDGEMENTS
The authors wish to thank Paul Stoffa, John Shepherd and Keith McCamy for many constructive criticisms. Funds for this work were pro¬ vided through the Climate Dynamics Program of the National Science Founda¬ tion, grant no. ATM-77-19217.
19
FIGURE CAPTIONS
Figure 1.
Figure 2.
Figure 3.
Figure 4.
Figure 5.
Figure 6 .
Figure 7.
The gain of the frequency response function u (f) for various
P
values of the parameter p over the range .001^f^.5.
The base ten logarithm of the negative of the frequency response function u ' (f) for various values of the parameter p over the
.c'
range . 001^f<.5. The dashed line passes through the maximi (circles) on these curves.
The root z++ of the characteristic polynomial r (z) . Values of P P
the parameter p are shown along the curve .
The central values vp(0) and v^ 1 (0) of the smooth parts of the impulse response functions u^ and u ' .
The smooth part v^(t) of the impulse response function u^ for three values of p most useful for standardizing forest interior tree-ring series.
The smooth part v ' (t) of the impulse response function u^ ' for three values of p most useful for standardizing forest interior tree-ring series.
Four examples of the smoothing spline applied to forest interior ring-width series. Each spline was computed for log p = -4.0. The solud curve is the spline fit and the dashed curve is the orthogonal polynomial fit from program INDXA (see ref. 4) . The four series are from the same site in southeastern New York and the lower two plots are curves from the same tree.
20
1.0
0.9
0.8
0.7 0.6
0.5' 0.4 0.3
02 0.1
0.0
-3.0 -2.0 -1.0 0.0
Q.
LOG f
Figure 1
2.5
LOG f
Figure 2
LOG [- up '( f )]
++
dZ
ro
Figure
LOG
-25
RINGWIDTHS (mm) RINGWIDTHS (mm)
1.5-
A
E
S
cn 1-0-
X
I-
9
5
o
z
* 0.5-
0.0
r— I— t— — t— i - 1 - r— t— — i - r - r - v. . , - r
1850 1900 1950
I 1 "—'-I 1 ■ " I ■ ' I "I1 ■ — T ' I - 1 I I . n - 1 ' 1 — 1111
1800 1850 1900 1950
— , - i 1 i - 1 i - 1 - 1 - r ■ 11 - 1 - 1 - 1 - 1 - 1 - 1 - - - 1 - 1 - 1 - -
1750 1800 1850 1900 1950
YEARS
Figure 7
REFERENCES
1. Churchill, R.V. : Introduction to complex variables and applications.
New York: McGraw-Hill Book Company, Inc., 1948.
2. Cook, E.R. and Peters, K.: The smoothing spline: a new approach to
standardizing forest interior ring-width series for dendroclimatic studies. (submitted) .
3. Craven, P.G.: Smoothing noisy data with spline functions, Numer. Math.
31, 377-403 (1979) .
4. Demmler, A., Reinsch, C.: Oscillation matrices with spline smoothing,
Numer. Math. 24, 375-382 (1975).
5. Fritts, H.C., Mosimann, J.E. and Bottorff, C.P.: A revised computer
program for standardizing tree-ring series. Tree-Ring Bulletin 29, 15-20 (1969).
6. Horowitz, L.L. : The effects of spline interpolation on power spectral
density. IEEE Transactions on Acoustics, Speech, and Signal Pro¬ cessing, vol. AAP-22 , number 1, 22-26 (February 1974) .
7. Reinsch, C.H. : Smoothing by spline functions, Numer. Math. 10, 177-183
(1967) .
8. Reinsch, C.H„: Smoothing by spline functions, II. Numer. Math. 16,
451-454 (1971) .
9. Wahba, G. : Smoothing noisy data with spline functions, Numer. Math. 24,
383-393 (1975) .
10. Wahba, G. , Wold, S.: A completely automatic French curve: Fitting
spline functions by cross-validation, Comm. Stat. 4, 1-17 (1975).
11. Wahba, G. , Wold, S.: Periodic splines for spectral density estimation:
The use of cross-validation for determining the degree of smoothing. Comm. Stat. 4, 125-141 (1975).
12. Wold, S.: Spline functions in data analysis, Technometrics 16, 1-11 (1974) .
21
LAMONT DOHERTY GEOLOGICAL OBSERVATORY
OF COLUMBIA UNIVERSITY